Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-12-20 03:14:04 

0001 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0002 #include <Math/VectorUtil.h>
0003 #include "DataFormats/Math/interface/deltaR.h"
0004 #include "DataFormats/HGCalReco/interface/Common.h"
0005 #include "DataFormats/GeometrySurface/interface/BoundDisk.h"
0006 #include "DataFormats/HGCalReco/interface/Trackster.h"
0007 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h"
0008 #include "RecoParticleFlow/PFProducer/interface/PFMuonAlgo.h"
0009 #include "RecoHGCal/TICL/interface/TracksterLinkingAlgoBase.h"
0010 #include "RecoHGCal/TICL/plugins/TracksterLinkingbySkeletons.h"
0011 #include "TICLGraph.h"
0012 
0013 namespace {
0014   bool isRoundTrackster(std::array<ticl::Vector, 3> skeleton) { return (skeleton[0].Z() == skeleton[2].Z()); }
0015 
0016   bool isGoodTrackster(const ticl::Trackster &trackster,
0017                        const std::array<ticl::Vector, 3> &skeleton,
0018                        const unsigned int min_num_lcs,
0019                        const float min_trackster_energy,
0020                        const float pca_quality_th) {
0021     bool isGood = false;
0022 
0023     if (isRoundTrackster(skeleton) or trackster.vertices().size() < min_num_lcs) {
0024       if (trackster.raw_energy() > min_trackster_energy) {
0025         auto const &eigenvalues = trackster.eigenvalues();
0026         auto const sum = std::accumulate(std::begin(eigenvalues), std::end(eigenvalues), 0.f);
0027         float pcaQuality = sum > 0.f ? eigenvalues[0] / sum : 0.f;
0028         if (pcaQuality > pca_quality_th) {
0029           isGood = true;
0030         }
0031       }
0032     } else {
0033       auto const &eigenvalues = trackster.eigenvalues();
0034       auto const sum = std::accumulate(std::begin(eigenvalues), std::end(eigenvalues), 0.f);
0035       float pcaQuality = sum > 0.f ? eigenvalues[0] / sum : 0.f;
0036       if (pcaQuality > pca_quality_th) {
0037         isGood = true;
0038       }
0039     }
0040     return isGood;
0041   }
0042 
0043   //distance between skeletons
0044   inline float projective_distance(const ticl::Vector &point1, const ticl::Vector &point2) {
0045     // squared projective distance
0046     float r1 = std::sqrt(point1.x() * point1.x() + point1.y() * point1.y());
0047     float r2_at_z1 =
0048         std::sqrt(point2.x() * point2.x() + point2.y() * point2.y()) * std::abs(point1.z()) / std::abs(point2.z());
0049     float delta_phi = reco::deltaPhi(point1.Phi(), point2.Phi());
0050     float projective_distance = (r1 - r2_at_z1) * (r1 - r2_at_z1) + r2_at_z1 * r2_at_z1 * delta_phi * delta_phi;
0051     LogDebug("TracksterLinkingbySkeletons") << "Computing distance between point : " << point1 << " And point "
0052                                             << point2 << " Distance " << projective_distance << std::endl;
0053     return projective_distance;
0054   }
0055 }  // namespace
0056 
0057 using namespace ticl;
0058 
0059 TracksterLinkingbySkeletons::TracksterLinkingbySkeletons(const edm::ParameterSet &conf,
0060                                                          edm::ConsumesCollector iC,
0061                                                          cms::Ort::ONNXRuntime const *onnxRuntime)
0062     : TracksterLinkingAlgoBase(conf, iC),
0063       lower_boundary_(conf.getParameter<std::vector<double>>("lower_boundary")),
0064       upper_boundary_(conf.getParameter<std::vector<double>>("upper_boundary")),
0065       upper_distance_projective_sqr_(conf.getParameter<std::vector<double>>("upper_distance_projective_sqr")),
0066       lower_distance_projective_sqr_(conf.getParameter<std::vector<double>>("lower_distance_projective_sqr")),
0067       min_distance_z_(conf.getParameter<std::vector<double>>("min_distance_z")),
0068       upper_distance_projective_sqr_closest_points_(
0069           conf.getParameter<std::vector<double>>("upper_distance_projective_sqr_closest_points")),
0070       lower_distance_projective_sqr_closest_points_(
0071           conf.getParameter<std::vector<double>>("lower_distance_projective_sqr_closest_points")),
0072       max_z_distance_closest_points_(conf.getParameter<std::vector<double>>("max_z_distance_closest_points")),
0073       cylinder_radius_sqr_(conf.getParameter<std::vector<double>>("cylinder_radius_sqr")),
0074       cylinder_radius_sqr_split_(conf.getParameter<double>("cylinder_radius_sqr_split")),
0075       proj_distance_split_(conf.getParameter<double>("proj_distance_split")),
0076       timing_quality_threshold_(conf.getParameter<double>("track_time_quality_threshold")),
0077       min_trackster_energy_(conf.getParameter<double>("min_trackster_energy")),
0078       pca_quality_th_(conf.getParameter<double>("pca_quality_th")),
0079       dot_prod_th_(conf.getParameter<double>("dot_prod_th")),
0080       deltaRxy_(conf.getParameter<double>("deltaRxy")),
0081       min_num_lcs_(conf.getParameter<unsigned int>("min_num_lcs"))
0082 
0083 {}
0084 
0085 void TracksterLinkingbySkeletons::buildLayers() {
0086   // build disks at HGCal front & EM-Had interface for track propagation
0087 
0088   float zVal = hgcons_->waferZ(1, true);
0089   std::pair<float, float> rMinMax = hgcons_->rangeR(zVal, true);
0090 
0091   float zVal_interface = rhtools_.getPositionLayer(rhtools_.lastLayerEE()).z();
0092   std::pair<float, float> rMinMax_interface = hgcons_->rangeR(zVal_interface, true);
0093 
0094   for (int iSide = 0; iSide < 2; ++iSide) {
0095     float zSide = (iSide == 0) ? (-1. * zVal) : zVal;
0096     firstDisk_[iSide] =
0097         std::make_unique<GeomDet>(Disk::build(Disk::PositionType(0, 0, zSide),
0098                                               Disk::RotationType(),
0099                                               SimpleDiskBounds(rMinMax.first, rMinMax.second, zSide - 0.5, zSide + 0.5))
0100                                       .get());
0101 
0102     zSide = (iSide == 0) ? (-1. * zVal_interface) : zVal_interface;
0103     interfaceDisk_[iSide] = std::make_unique<GeomDet>(
0104         Disk::build(Disk::PositionType(0, 0, zSide),
0105                     Disk::RotationType(),
0106                     SimpleDiskBounds(rMinMax_interface.first, rMinMax_interface.second, zSide - 0.5, zSide + 0.5))
0107             .get());
0108   }
0109 }
0110 
0111 void TracksterLinkingbySkeletons::initialize(const HGCalDDDConstants *hgcons,
0112                                              const hgcal::RecHitTools rhtools,
0113                                              const edm::ESHandle<MagneticField> bfieldH,
0114                                              const edm::ESHandle<Propagator> propH) {
0115   hgcons_ = hgcons;
0116   rhtools_ = rhtools;
0117   buildLayers();
0118 
0119   bfield_ = bfieldH;
0120   propagator_ = propH;
0121 
0122   //define LUT for eta windows
0123   // eta windows obtained with a deltaR of 4cm at z = 400 cm
0124   float etaStep = (TileConstants::maxEta - TileConstants::minEta) / TileConstants::nEtaBins;
0125   float expNeg2DeltaRxy = deltaRxy_ * std::exp(-2.f);
0126 
0127   for (int i = 0; i < TileConstants::nEtaBins; ++i) {
0128     float eta = TileConstants::minEta + i * etaStep;
0129 
0130     float expNegEta = std::exp(-eta);
0131     float R = z_surface * 2.f * expNegEta / (1.f - expNeg2DeltaRxy * expNegEta);
0132 
0133     eta_windows_[i] = std::abs(atan(deltaRxy_ / R));
0134   }
0135 }
0136 
0137 std::array<ticl::Vector, 3> TracksterLinkingbySkeletons::findSkeletonNodes(
0138     const ticl::Trackster &trackster,
0139     float lower_percentage,
0140     float upper_percentage,
0141     const std::vector<reco::CaloCluster> &layerClusters,
0142     const hgcal::RecHitTools &rhtools) {
0143   auto const &vertices = trackster.vertices();
0144   auto const trackster_raw_energy = trackster.raw_energy();
0145   // sort vertices by layerId
0146   std::array<ticl::Vector, 3> skeleton;
0147   if (trackster.vertices().size() < 3) {
0148     const auto &v = layerClusters[trackster.vertices()[0]];
0149     const Vector intersection(v.x(), v.y(), v.z());
0150     skeleton = {{intersection, intersection, intersection}};
0151     return skeleton;
0152   }
0153 
0154   std::vector<unsigned int> sortedVertices(vertices);
0155   std::sort(sortedVertices.begin(), sortedVertices.end(), [&layerClusters](unsigned int i, unsigned int j) {
0156     return std::abs(layerClusters[i].z()) < std::abs(layerClusters[j].z());
0157   });
0158 
0159   // now loop over sortedVertices and find the layerId that contains the lower_percentage of the energy
0160   // and the layerId that contains the upper_percentage of the energy
0161   float cumulativeEnergyFraction = 0.f;
0162   int innerLayerId = rhtools.getLayerWithOffset(layerClusters[sortedVertices[0]].hitsAndFractions()[0].first);
0163   float innerLayerZ = layerClusters[sortedVertices[0]].z();
0164   int outerLayerId = rhtools.getLayerWithOffset(layerClusters[sortedVertices.back()].hitsAndFractions()[0].first);
0165   float outerLayerZ = layerClusters[sortedVertices.back()].z();
0166   bool foundInnerLayer = false;
0167   bool foundOuterLayer = false;
0168   for (auto const &v : sortedVertices) {
0169     auto const &lc = layerClusters[v];
0170     auto const &n_lay = rhtools.getLayerWithOffset(lc.hitsAndFractions()[0].first);
0171     cumulativeEnergyFraction += lc.energy() / trackster_raw_energy;
0172     if (cumulativeEnergyFraction >= lower_percentage and not foundInnerLayer) {
0173       innerLayerId = n_lay;
0174       innerLayerZ = lc.z();
0175       foundInnerLayer = true;
0176     }
0177     if (cumulativeEnergyFraction >= upper_percentage and not foundOuterLayer) {
0178       outerLayerId = n_lay;
0179       outerLayerZ = lc.z();
0180       foundOuterLayer = true;
0181     }
0182   }
0183   int minimumDistanceInLayers = 4;
0184   if (outerLayerId - innerLayerId < minimumDistanceInLayers) {
0185     skeleton = {{trackster.barycenter(), trackster.barycenter(), trackster.barycenter()}};
0186   } else {
0187     auto intersectLineWithSurface = [](float surfaceZ, const Vector &origin, const Vector &direction) -> Vector {
0188       auto const t = (surfaceZ - origin.Z()) / direction.Z();
0189       auto const iX = t * direction.X() + origin.X();
0190       auto const iY = t * direction.Y() + origin.Y();
0191       auto const iZ = surfaceZ;
0192       const Vector intersection(iX, iY, iZ);
0193       return intersection;
0194     };
0195 
0196     auto const &t0_p1 = trackster.barycenter();
0197     auto const t0_p0 = intersectLineWithSurface(innerLayerZ, t0_p1, trackster.eigenvectors(0));
0198     auto const t0_p2 = intersectLineWithSurface(outerLayerZ, t0_p1, trackster.eigenvectors(0));
0199     skeleton = {{t0_p0, t0_p1, t0_p2}};
0200     std::sort(skeleton.begin(), skeleton.end(), [](Vector &v1, Vector &v2) { return v1.Z() < v2.Z(); });
0201   }
0202 
0203   return skeleton;
0204 }
0205 
0206 inline bool isInCylinder(const std::array<ticl::Vector, 3> &mySkeleton,
0207                          const std::array<ticl::Vector, 3> &otherSkeleton,
0208                          const float radius_sqr) {
0209   const auto &first = mySkeleton[0];
0210   const auto &last = mySkeleton[2];
0211   const auto &pointToCheck = otherSkeleton[0];
0212 
0213   const auto &cylAxis = last - first;
0214   const auto &vecToPoint = pointToCheck - first;
0215 
0216   auto axisNorm = cylAxis.Dot(cylAxis);
0217   auto projLength = vecToPoint.Dot(cylAxis) / axisNorm;
0218   bool isWithinLength = projLength >= 0 && projLength <= 1;
0219 
0220   if (!isWithinLength)
0221     return false;
0222 
0223   const auto &proj = cylAxis * projLength;
0224 
0225   const auto &pointOnAxis = first + proj;
0226 
0227   const auto &distance = pointToCheck - pointOnAxis;
0228   auto distance2 = distance.Dot(distance);
0229   LogDebug("TracksterLinkingbySkeletons") << "is within lenght " << distance2 << " radius " << radius_sqr << std::endl;
0230   bool isWithinRadius = distance2 <= radius_sqr;
0231 
0232   return isWithinRadius;
0233 }
0234 
0235 inline float computeParameter(float energy, float en_th_low, float cut1, float en_th_high, float cut2) {
0236   if (en_th_low == en_th_high) {  //protect if the thresholds are equal
0237     return (energy <= en_th_low) ? cut1 : cut2;
0238   }
0239   if (energy < en_th_low) {
0240     return cut1;
0241   } else if (energy >= en_th_low && energy <= en_th_high) {
0242     return ((cut2 - cut1) / (en_th_high - en_th_low)) * (energy - en_th_low) + cut1;
0243   } else {
0244     return cut2;
0245   }
0246 }
0247 
0248 // function to check if otherTrackster is a splitted component of myTrackster, meaning that otherTrackster is very close to myTrackster and seems to be generated by splitting of 2D layer clusters
0249 inline bool TracksterLinkingbySkeletons::isSplitComponent(const ticl::Trackster &myTrackster,
0250                                                           const ticl::Trackster &otherTrackster,
0251                                                           const std::array<ticl::Vector, 3> &mySkeleton,
0252                                                           const std::array<ticl::Vector, 3> &otherSkeleton,
0253                                                           float proj_distance) {
0254   //check if otherSKeleton z is within the z range of mySkeleton
0255   if (otherSkeleton[1].z() < mySkeleton[2].z() && otherSkeleton[1].z() > mySkeleton[0].z()) {
0256     if (proj_distance < proj_distance_split_) {
0257       return true;
0258     } else {
0259       // check if barycenter of otherTrackster is within the cirlce of 3cm of myTrackster barycenter
0260       //compute XY distance between barycenters
0261       float distance2 = (myTrackster.barycenter().x() - otherTrackster.barycenter().x()) *
0262                             (myTrackster.barycenter().x() - otherTrackster.barycenter().x()) +
0263                         (myTrackster.barycenter().y() - otherTrackster.barycenter().y()) *
0264                             (myTrackster.barycenter().y() - otherTrackster.barycenter().y());
0265       if (distance2 < cylinder_radius_sqr_split_) {
0266         return true;
0267       }
0268     }
0269   }
0270   return false;
0271 }
0272 
0273 bool TracksterLinkingbySkeletons::areCompatible(const ticl::Trackster &myTrackster,
0274                                                 const ticl::Trackster &otherTrackster,
0275                                                 const std::array<ticl::Vector, 3> &mySkeleton,
0276                                                 const std::array<ticl::Vector, 3> &otherSkeleton) {
0277   float zVal_interface = rhtools_.getPositionLayer(rhtools_.lastLayerEE()).z();
0278 
0279   if (!isGoodTrackster(myTrackster, mySkeleton, min_num_lcs_, min_trackster_energy_, pca_quality_th_)) {
0280     LogDebug("TracksterLinkingbySkeletons") << "Inner Trackster with energy " << myTrackster.raw_energy() << " Num LCs "
0281                                             << myTrackster.vertices().size() << " NOT GOOD " << std::endl;
0282     return false;
0283   }
0284 
0285   LogDebug("TracksterLinkingbySkeletons") << "Inner Trackster with energy " << myTrackster.raw_energy() << " Num LCs "
0286                                           << myTrackster.vertices().size() << " IS GOOD " << std::endl;
0287 
0288   float proj_distance = projective_distance(mySkeleton[1], otherSkeleton[1]);
0289   auto isEE = mySkeleton[1].z() <= zVal_interface ? 0 : 1;
0290   auto const max_distance_proj_sqr = computeParameter(myTrackster.raw_energy(),
0291                                                       lower_boundary_[isEE],
0292                                                       lower_distance_projective_sqr_[isEE],
0293                                                       upper_boundary_[isEE],
0294                                                       upper_distance_projective_sqr_[isEE]);
0295   bool areAlignedInProjectiveSpace = proj_distance < max_distance_proj_sqr;
0296 
0297   LogDebug("TracksterLinkingbySkeletons")
0298       << "\t Trying to compare with outer Trackster with energy " << otherTrackster.raw_energy() << " Num LCS "
0299       << otherTrackster.vertices().size() << " Projective distance " << proj_distance << " areAlignedProjective "
0300       << areAlignedInProjectiveSpace << " TH " << max_distance_proj_sqr << std::endl;
0301 
0302   if (isGoodTrackster(otherTrackster, otherSkeleton, min_num_lcs_, min_trackster_energy_, pca_quality_th_)) {
0303     if (areAlignedInProjectiveSpace) {
0304       LogDebug("TracksterLinkingbySkeletons") << "\t\t Linked! " << std::endl;
0305       return true;
0306     } else {
0307       //if the tracksters are not aligned in Projective distance, check if otherTrackster is within the cylinder of 3cm radius
0308       //this is used to recover LC splittings
0309       if (isSplitComponent(myTrackster, otherTrackster, mySkeleton, otherSkeleton, proj_distance)) {
0310         LogDebug("TracksterLinkingbySkeletons") << "\t\t Linked! Splitted components!" << std::endl;
0311         return true;
0312       }
0313       //if is EE do not try to link more, PU occupancy is too high in this region
0314       if (isEE) {
0315         return false;
0316       }
0317       //if instead we are in the CE-H part of the detector, we can try to link more
0318       // we measure the distance between the two closest nodes in the two skeletons
0319       return checkClosestPoints(myTrackster, otherTrackster, mySkeleton, otherSkeleton, isEE);
0320     }
0321   } else {
0322     if (otherTrackster.vertices().size() >= 3) {
0323       if (areAlignedInProjectiveSpace) {
0324         LogDebug("TracksterLinkingbySkeletons") << "\t\t Linked! " << std::endl;
0325         return true;
0326       } else {
0327         LogDebug("TracksterLinkingbySkeletons")
0328             << "\t Not aligned in projective space, check distance between closest points in the two skeletons "
0329             << std::endl;
0330         if (checkClosestPoints(myTrackster, otherTrackster, mySkeleton, otherSkeleton, isEE)) {
0331           return true;
0332         } else {
0333           return checkCylinderAlignment(mySkeleton, otherSkeleton, isEE);
0334         }
0335       }
0336     } else {
0337       return checkCylinderAlignment(mySkeleton, otherSkeleton, isEE);
0338     }
0339   }
0340 }
0341 
0342 bool TracksterLinkingbySkeletons::checkCylinderAlignment(const std::array<ticl::Vector, 3> &mySkeleton,
0343                                                          const std::array<ticl::Vector, 3> &otherSkeleton,
0344                                                          int isEE) {
0345   bool isInCyl = isInCylinder(mySkeleton, otherSkeleton, cylinder_radius_sqr_[isEE]);
0346   if (isInCyl) {
0347     LogDebug("TracksterLinkingbySkeletons") << "Two Points are in Cylinder  " << isInCyl << " Linked! " << std::endl;
0348   }
0349   return isInCyl;
0350 }
0351 
0352 bool TracksterLinkingbySkeletons::checkSkeletonAlignment(const ticl::Trackster &myTrackster,
0353                                                          const ticl::Trackster &otherTrackster) {
0354   auto dotProdSkeletons = myTrackster.eigenvectors(0).Dot(otherTrackster.eigenvectors(0));
0355   bool alignedSkeletons = dotProdSkeletons > dot_prod_th_;
0356 
0357   LogDebug("TracksterLinkingbySkeletons")
0358       << "\t Outer Trackster is Good, checking for skeleton alignment " << alignedSkeletons << " dotProd "
0359       << dotProdSkeletons << " Threshold " << dot_prod_th_ << std::endl;
0360 
0361   if (alignedSkeletons) {
0362     LogDebug("TracksterLinkingbySkeletons") << "\t\t Linked! " << std::endl;
0363   }
0364 
0365   return alignedSkeletons;
0366 }
0367 
0368 bool TracksterLinkingbySkeletons::checkClosestPoints(const ticl::Trackster &myTrackster,
0369                                                      const ticl::Trackster &otherTrackster,
0370                                                      const std::array<ticl::Vector, 3> &mySkeleton,
0371                                                      const std::array<ticl::Vector, 3> &otherSkeleton,
0372                                                      int isEE) {
0373   int myClosestPoint = -1;
0374   int otherClosestPoint = -1;
0375   float minDistance_z = std::numeric_limits<float>::max();
0376 
0377   for (int i = 1; i < 3; i++) {
0378     for (int j = 0; j < 3; j++) {
0379       float dist_z = std::abs(mySkeleton[i].Z() - otherSkeleton[j].Z());
0380       if (dist_z < minDistance_z) {
0381         myClosestPoint = i;
0382         otherClosestPoint = j;
0383         minDistance_z = dist_z;
0384       }
0385     }
0386   }
0387 
0388   float d = projective_distance(mySkeleton[myClosestPoint], otherSkeleton[otherClosestPoint]);
0389   auto const max_distance_proj_sqr_closest = computeParameter(myTrackster.raw_energy(),
0390                                                               lower_boundary_[isEE],
0391                                                               lower_distance_projective_sqr_closest_points_[isEE],
0392                                                               upper_boundary_[isEE],
0393                                                               upper_distance_projective_sqr_closest_points_[isEE]);
0394 
0395   LogDebug("TracksterLinkingbySkeletons")
0396       << "\t\t Distance between closest points " << d << " TH " << 10.f << " Z Distance " << minDistance_z << " TH "
0397       << max_distance_proj_sqr_closest << std::endl;
0398 
0399   if (d < max_distance_proj_sqr_closest && minDistance_z < max_z_distance_closest_points_[isEE]) {
0400     LogDebug("TracksterLinkingbySkeletons") << "\t\t\t Linked! " << d << std::endl;
0401     return true;
0402   }
0403 
0404   LogDebug("TracksterLinkingbySkeletons") << "Distance between closest point " << d << " Distance in z "
0405                                           << max_z_distance_closest_points_[isEE] << std::endl;
0406 
0407   return checkCylinderAlignment(mySkeleton, otherSkeleton, isEE);
0408 }
0409 
0410 void TracksterLinkingbySkeletons::linkTracksters(
0411     const Inputs &input,
0412     std::vector<Trackster> &resultTracksters,
0413     std::vector<std::vector<unsigned int>> &linkedResultTracksters,
0414     std::vector<std::vector<unsigned int>> &linkedTracksterIdToInputTracksterId) {
0415   const auto &tracksters = input.tracksters;
0416   const auto &layerClusters = input.layerClusters;
0417 
0418   // sort tracksters by energy
0419   std::vector<unsigned int> sortedTracksters(tracksters.size());
0420   std::iota(sortedTracksters.begin(), sortedTracksters.end(), 0);
0421   std::sort(sortedTracksters.begin(), sortedTracksters.end(), [&tracksters](unsigned int i, unsigned int j) {
0422     return tracksters[i].raw_energy() > tracksters[j].raw_energy();
0423   });
0424   // fill tiles for trackster linking
0425   // tile 0 for negative eta
0426   // tile 1 for positive eta
0427   std::array<TICLLayerTile, 2> tracksterTile;
0428   // loop over tracksters sorted by energy and calculate skeletons
0429   // fill tiles for trackster linking
0430   std::vector<std::array<ticl::Vector, 3>> skeletons(tracksters.size());
0431   for (auto const t_idx : sortedTracksters) {
0432     const auto &trackster = tracksters[t_idx];
0433     skeletons[t_idx] = findSkeletonNodes(tracksters[t_idx], 0.1, 0.9, layerClusters, rhtools_);
0434     tracksterTile[trackster.barycenter().eta() > 0.f].fill(
0435         trackster.barycenter().eta(), trackster.barycenter().phi(), t_idx);
0436   }
0437   std::vector<int> maskReceivedLink(tracksters.size(), 1);
0438   std::vector<int> isRootTracksters(tracksters.size(), 1);
0439 
0440   std::vector<ticl::Node> allNodes;
0441   for (size_t it = 0; it < tracksters.size(); ++it) {
0442     allNodes.emplace_back(it);
0443   }
0444 
0445   // loop over tracksters sorted by energy and link them
0446   for (auto const &t_idx : sortedTracksters) {
0447     auto const &trackster = tracksters[t_idx];
0448     auto const &skeleton = skeletons[t_idx];
0449 
0450     auto const &bary = trackster.barycenter();
0451     int tileIndex = bary.eta() > 0.f;
0452     const auto &tiles = tracksterTile[tileIndex];
0453     auto const window = eta_windows_[tiles.etaBin(bary.eta())];
0454     float eta_min = std::max(abs(bary.eta()) - window, TileConstants::minEta);
0455     float eta_max = std::min(abs(bary.eta()) + window, TileConstants::maxEta);
0456     std::array<int, 4> search_box = tiles.searchBoxEtaPhi(eta_min, eta_max, bary.phi() - window, bary.phi() + window);
0457 
0458     if (search_box[2] > search_box[3]) {
0459       search_box[3] += TileConstants::nPhiBins;
0460     }
0461 
0462     for (int eta_i = search_box[0]; eta_i <= search_box[1]; ++eta_i) {
0463       for (int phi_i = search_box[2]; phi_i <= search_box[3]; ++phi_i) {
0464         auto &neighbours = tiles[tiles.globalBin(eta_i, (phi_i % TileConstants::nPhiBins))];
0465         for (auto n : neighbours) {
0466           if (t_idx == n)
0467             continue;
0468 
0469           auto const &tracksterOut = tracksters[n];
0470           auto const &skeletonOut = skeletons[n];
0471           auto const deltaphi = reco::deltaPhi(trackster.barycenter().phi(), tracksterOut.barycenter().phi());
0472           if (abs(trackster.barycenter().eta() - tracksterOut.barycenter().eta()) <= window && deltaphi <= window) {
0473             bool isInGood = isGoodTrackster(trackster, skeleton, min_num_lcs_, min_trackster_energy_, pca_quality_th_);
0474             bool isOutGood =
0475                 isGoodTrackster(tracksterOut, skeletonOut, min_num_lcs_, min_trackster_energy_, pca_quality_th_);
0476             if (isInGood) {
0477               LogDebug("TracksterLinkingbySkeletons")
0478                   << "Trying to Link Trackster " << t_idx << " With Trackster " << n << std::endl;
0479               if (areCompatible(trackster, tracksters[n], skeleton, skeletonOut)) {
0480                 LogDebug("TracksterLinkingbySkeletons")
0481                     << "\t==== LINK: Trackster " << t_idx << " Linked with Trackster " << n << std::endl;
0482                 //    maskReceivedLink[n] = 0;
0483                 if (isOutGood) {
0484                   if (abs(skeleton[0].z()) < abs(skeletonOut[0].z())) {
0485                     LogDebug("TracksterLinkingbySkeletons") << "Trackster in energy " << trackster.raw_energy()
0486                                                             << " Out is good " << tracksterOut.raw_energy() << " Sk In "
0487                                                             << skeleton[0] << " Sk out " << skeletonOut[0] << std::endl;
0488                     LogDebug("TracksterLinkingbySkeletons") << "\t " << t_idx << " --> " << n << std::endl;
0489                     allNodes[t_idx].addOuterNeighbour(n);
0490                     allNodes[n].addInnerNeighbour(t_idx);
0491                     isRootTracksters[n] = 0;
0492                   } else if (abs(skeleton[0].z()) > abs(skeletonOut[0].z())) {
0493                     LogDebug("TracksterLinkingbySkeletons") << "Trackster in energy " << trackster.raw_energy()
0494                                                             << " Out is good " << tracksterOut.raw_energy() << " Sk In "
0495                                                             << skeleton[0] << " Sk out " << skeletonOut[0] << std::endl;
0496                     LogDebug("TracksterLinkingbySkeletons") << "\t " << n << " --> " << t_idx << std::endl;
0497                     allNodes[n].addOuterNeighbour(t_idx);
0498                     allNodes[t_idx].addInnerNeighbour(n);
0499                     isRootTracksters[t_idx] = 0;
0500                   } else {
0501                     if (trackster.raw_energy() >= tracksterOut.raw_energy()) {
0502                       LogDebug("TracksterLinkingbySkeletons")
0503                           << "Trackster in energy " << trackster.raw_energy() << " Out is good "
0504                           << tracksterOut.raw_energy() << " Sk In " << skeleton[0] << " Sk out " << skeletonOut[0]
0505                           << std::endl;
0506                       LogDebug("TracksterLinkingbySkeletons") << "\t " << t_idx << " --> " << n << std::endl;
0507                       allNodes[t_idx].addOuterNeighbour(n);
0508                       allNodes[n].addInnerNeighbour(t_idx);
0509                       isRootTracksters[n] = 0;
0510                     } else {
0511                       LogDebug("TracksterLinkingbySkeletons")
0512                           << "Trackster in energy " << trackster.raw_energy() << " Out is good "
0513                           << tracksterOut.raw_energy() << " Sk In " << skeleton[0] << " Sk out " << skeletonOut[0]
0514                           << std::endl;
0515                       LogDebug("TracksterLinkingbySkeletons") << "\t " << n << " --> " << t_idx << std::endl;
0516                       allNodes[n].addOuterNeighbour(t_idx);
0517                       allNodes[t_idx].addInnerNeighbour(n);
0518                       isRootTracksters[t_idx] = 0;
0519                     }
0520                   }
0521                 } else {
0522                   LogDebug("TracksterLinkingbySkeletons")
0523                       << "Trackster in energy " << trackster.raw_energy() << " Out is NOT good "
0524                       << tracksterOut.raw_energy() << " Sk In " << skeleton[0] << " Sk out " << skeletonOut[0]
0525                       << std::endl;
0526                   LogDebug("TracksterLinkingbySkeletons") << "\t " << t_idx << " --> " << n << std::endl;
0527                   allNodes[t_idx].addOuterNeighbour(n);
0528                   allNodes[n].addInnerNeighbour(t_idx);
0529                   isRootTracksters[n] = 0;
0530                 }
0531               }
0532             }
0533           }
0534         }
0535       }
0536     }
0537   }
0538 
0539   TICLGraph graph(allNodes);
0540   auto sortedRootNodes = graph.getRootNodes();
0541   std::sort(sortedRootNodes.begin(), sortedRootNodes.end(), [&tracksters](const ticl::Node &n1, const ticl::Node &n2) {
0542     unsigned int n1Id = n1.getId();
0543     unsigned int n2Id = n2.getId();
0544     return tracksters[n1Id].raw_energy() > tracksters[n2Id].raw_energy();
0545   });
0546 
0547   int ic = 0;
0548   auto const &components = graph.findSubComponents(sortedRootNodes);
0549   linkedTracksterIdToInputTracksterId.resize(components.size());
0550   for (auto const &comp : components) {
0551     LogDebug("TracksterLinkingbySkeletons") << "Component " << ic << " Node: ";
0552     std::vector<unsigned int> linkedTracksters;
0553     Trackster outTrackster;
0554     if (comp.size() == 1) {
0555       if (input.tracksters[comp[0]].vertices().size() <= 3 && input.tracksters[comp[0]].raw_energy() < 5.f) {
0556         continue;
0557       }
0558     }
0559     for (auto const &node : comp) {
0560       LogDebug("TracksterLinkingbySkeletons") << node << " ";
0561       linkedTracksterIdToInputTracksterId[ic].push_back(node);
0562       outTrackster.mergeTracksters(input.tracksters[node]);
0563     }
0564     linkedTracksters.push_back(resultTracksters.size());
0565     LogDebug("TracksterLinkingbySkeletons") << "\nOut Trackster " << outTrackster.raw_energy() << std::endl;
0566     resultTracksters.push_back(outTrackster);
0567     linkedResultTracksters.push_back(linkedTracksters);
0568     LogDebug("TracksterLinkingbySkeletons") << "\n";
0569     ++ic;
0570   }
0571   LogDebug("TracksterLinkingbySkeletons") << "\n";
0572 
0573 }  // linkTracksters
[ 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!