Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2024-09-26 05:07:08

0001 // Author: Marco Rovere - marco.rovere@cern.ch
0002 // Date: 04/2021
0003 #include <algorithm>
0004 #include <set>
0005 #include <vector>
0006 
0007 #include "DataFormats/Math/interface/deltaR.h"
0008 #include "FWCore/Framework/interface/Event.h"
0009 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0010 #include "FWCore/Utilities/interface/Exception.h"
0011 #include "PatternRecognitionbyCLUE3D.h"
0012 
0013 #include "TrackstersPCA.h"
0014 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
0015 #include "Geometry/Records/interface/CaloGeometryRecord.h"
0016 #include "FWCore/Framework/interface/EventSetup.h"
0017 
0018 using namespace ticl;
0019 
0020 template <typename TILES>
0021 PatternRecognitionbyCLUE3D<TILES>::PatternRecognitionbyCLUE3D(const edm::ParameterSet &conf, edm::ConsumesCollector iC)
0022     : PatternRecognitionAlgoBaseT<TILES>(conf, iC),
0023       caloGeomToken_(iC.esConsumes<CaloGeometry, CaloGeometryRecord>()),
0024       criticalDensity_(conf.getParameter<std::vector<double>>("criticalDensity")),
0025       criticalSelfDensity_(conf.getParameter<std::vector<double>>("criticalSelfDensity")),
0026       densitySiblingLayers_(conf.getParameter<std::vector<int>>("densitySiblingLayers")),
0027       densityEtaPhiDistanceSqr_(conf.getParameter<std::vector<double>>("densityEtaPhiDistanceSqr")),
0028       densityXYDistanceSqr_(conf.getParameter<std::vector<double>>("densityXYDistanceSqr")),
0029       kernelDensityFactor_(conf.getParameter<std::vector<double>>("kernelDensityFactor")),
0030       densityOnSameLayer_(conf.getParameter<bool>("densityOnSameLayer")),
0031       nearestHigherOnSameLayer_(conf.getParameter<bool>("nearestHigherOnSameLayer")),
0032       useAbsoluteProjectiveScale_(conf.getParameter<bool>("useAbsoluteProjectiveScale")),
0033       useClusterDimensionXY_(conf.getParameter<bool>("useClusterDimensionXY")),
0034       rescaleDensityByZ_(conf.getParameter<bool>("rescaleDensityByZ")),
0035       criticalEtaPhiDistance_(conf.getParameter<std::vector<double>>("criticalEtaPhiDistance")),
0036       criticalXYDistance_(conf.getParameter<std::vector<double>>("criticalXYDistance")),
0037       criticalZDistanceLyr_(conf.getParameter<std::vector<int>>("criticalZDistanceLyr")),
0038       outlierMultiplier_(conf.getParameter<std::vector<double>>("outlierMultiplier")),
0039       minNumLayerCluster_(conf.getParameter<std::vector<int>>("minNumLayerCluster")),
0040       doPidCut_(conf.getParameter<bool>("doPidCut")),
0041       cutHadProb_(conf.getParameter<double>("cutHadProb")),
0042       computeLocalTime_(conf.getParameter<bool>("computeLocalTime")),
0043       usePCACleaning_(conf.getParameter<bool>("usePCACleaning")){};
0044 template <typename TILES>
0045 void PatternRecognitionbyCLUE3D<TILES>::dumpTiles(const TILES &tiles) const {
0046   constexpr int nEtaBin = TILES::constants_type_t::nEtaBins;
0047   constexpr int nPhiBin = TILES::constants_type_t::nPhiBins;
0048   auto lastLayerPerSide = static_cast<int>(rhtools_.lastLayer(false));
0049   int maxLayer = 2 * lastLayerPerSide - 1;
0050   for (int layer = 0; layer <= maxLayer; layer++) {
0051     for (int ieta = 0; ieta < nEtaBin; ieta++) {
0052       auto offset = ieta * nPhiBin;
0053       for (int phi = 0; phi < nPhiBin; phi++) {
0054         int iphi = ((phi % nPhiBin + nPhiBin) % nPhiBin);
0055         if (!tiles[layer][offset + iphi].empty()) {
0056           if (this->algo_verbosity_ > VerbosityLevel::Advanced) {
0057             edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "Layer: " << layer << " ieta: " << ieta << " phi: " << phi
0058                                                            << " " << tiles[layer][offset + iphi].size();
0059           }
0060         }
0061       }
0062     }
0063   }
0064 }
0065 
0066 template <typename TILES>
0067 void PatternRecognitionbyCLUE3D<TILES>::dumpTracksters(const std::vector<std::pair<int, int>> &layerIdx2layerandSoa,
0068                                                        const int eventNumber,
0069                                                        const std::vector<Trackster> &tracksters) const {
0070   if (PatternRecognitionAlgoBaseT<TILES>::algo_verbosity_ > VerbosityLevel::Advanced) {
0071     edm::LogVerbatim("PatternRecognitionbyCLUE3D")
0072         << "[evt, tracksterId, cells, prob_photon, prob_ele, prob_chad, prob_nhad, layer_i, x_i, y_i, eta_i, phi_i, "
0073            "energy_i, radius_i, rho_i, z_extension, delta_tr, delta_lyr, isSeed_i";
0074   }
0075 
0076   int num = 0;
0077   const std::string sep(", ");
0078   for (auto const &t : tracksters) {
0079     for (auto v : t.vertices()) {
0080       auto [lyrIdx, soaIdx] = layerIdx2layerandSoa[v];
0081       auto const &thisLayer = clusters_[lyrIdx];
0082       if (PatternRecognitionAlgoBaseT<TILES>::algo_verbosity_ > VerbosityLevel::Advanced) {
0083         edm::LogVerbatim("PatternRecognitionbyCLUE3D_NTP")
0084             << std::setw(4) << eventNumber << sep << std::setw(4) << num << sep << std::setw(4) << t.vertices().size()
0085             << sep << std::setw(8) << t.id_probability(ticl::Trackster::ParticleType::photon) << sep << std::setw(8)
0086             << t.id_probability(ticl::Trackster::ParticleType::electron) << sep << std::setw(8)
0087             << t.id_probability(ticl::Trackster::ParticleType::charged_hadron) << sep << std::setw(8)
0088             << t.id_probability(ticl::Trackster::ParticleType::neutral_hadron) << sep << std::setw(4) << lyrIdx << sep
0089             << std::setw(10) << thisLayer.x[soaIdx] << sep << std::setw(10) << thisLayer.y[soaIdx] << sep
0090             << std::setw(10) << thisLayer.eta[soaIdx] << sep << std::setw(10) << thisLayer.phi[soaIdx] << sep
0091             << std::setw(10) << thisLayer.energy[soaIdx] << sep << std::setw(10) << thisLayer.radius[soaIdx] << sep
0092             << std::setw(10) << thisLayer.rho[soaIdx] << sep << std::setw(10) << thisLayer.z_extension[soaIdx] << sep
0093             << std::setw(10) << thisLayer.delta[soaIdx].first << sep << std::setw(10) << thisLayer.delta[soaIdx].second
0094             << sep << std::setw(4) << thisLayer.isSeed[soaIdx];
0095       }
0096     }
0097     num++;
0098   }
0099 }
0100 
0101 template <typename TILES>
0102 void PatternRecognitionbyCLUE3D<TILES>::dumpClusters(const TILES &tiles,
0103                                                      const std::vector<std::pair<int, int>> &layerIdx2layerandSoa,
0104                                                      const int eventNumber) const {
0105   if (PatternRecognitionAlgoBaseT<TILES>::algo_verbosity_ > VerbosityLevel::Advanced) {
0106     edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "[evt, lyr, Seed,      x,       y,       z, r/|z|,   eta,   phi, "
0107                                                       "etab,  phib, cells, enrgy, e/rho,   rho,   z_ext, "
0108                                                       "   dlt_tr,   dlt_lyr, "
0109                                                       " nestHL, nestHSoaIdx, radius, clIdx, lClOrigIdx, SOAidx";
0110   }
0111 
0112   for (unsigned int layer = 0; layer < clusters_.size(); layer++) {
0113     auto const &thisLayer = clusters_[layer];
0114     int num = 0;
0115     for (auto v : thisLayer.x) {
0116       if (PatternRecognitionAlgoBaseT<TILES>::algo_verbosity_ > VerbosityLevel::Advanced) {
0117         edm::LogVerbatim("PatternRecognitionbyCLUE3D")
0118             << std::setw(4) << eventNumber << ", " << std::setw(3) << layer << ", " << std::setw(4)
0119             << thisLayer.isSeed[num] << ", " << std::setprecision(3) << std::fixed << v << ", " << thisLayer.y[num]
0120             << ", " << thisLayer.z[num] << ", " << thisLayer.r_over_absz[num] << ", " << thisLayer.eta[num] << ", "
0121             << thisLayer.phi[num] << ", " << std::setw(5) << tiles[layer].etaBin(thisLayer.eta[num]) << ", "
0122             << std::setw(5) << tiles[layer].phiBin(thisLayer.phi[num]) << ", " << std::setw(4) << thisLayer.cells[num]
0123             << ", " << std::setprecision(3) << thisLayer.energy[num] << ", "
0124             << (thisLayer.energy[num] / thisLayer.rho[num]) << ", " << thisLayer.rho[num] << ", "
0125             << thisLayer.z_extension[num] << ", " << std::scientific << thisLayer.delta[num].first << ", "
0126             << std::setw(10) << thisLayer.delta[num].second << ", " << std::setw(5)
0127             << thisLayer.nearestHigher[num].first << ", " << std::setw(10) << thisLayer.nearestHigher[num].second
0128             << ", " << std::defaultfloat << std::setprecision(3) << thisLayer.radius[num] << ", " << std::setw(5)
0129             << thisLayer.clusterIndex[num] << ", " << std::setw(4) << thisLayer.layerClusterOriginalIdx[num] << ", "
0130             << std::setw(4) << num << ", ClusterInfo";
0131       }
0132       ++num;
0133     }
0134   }
0135   for (unsigned int lcIdx = 0; lcIdx < layerIdx2layerandSoa.size(); lcIdx++) {
0136     auto const &layerandSoa = layerIdx2layerandSoa[lcIdx];
0137     // Skip masked layer clusters
0138     if ((layerandSoa.first == -1) && (layerandSoa.second == -1))
0139       continue;
0140     if (PatternRecognitionAlgoBaseT<TILES>::algo_verbosity_ > VerbosityLevel::Advanced) {
0141       edm::LogVerbatim("PatternRecognitionbyCLUE3D")
0142           << "lcIdx: " << lcIdx << " on Layer: " << layerandSoa.first << " SOA: " << layerandSoa.second;
0143     }
0144   }
0145 }
0146 
0147 template <typename TILES>
0148 void PatternRecognitionbyCLUE3D<TILES>::makeTracksters(
0149     const typename PatternRecognitionAlgoBaseT<TILES>::Inputs &input,
0150     std::vector<Trackster> &result,
0151     std::unordered_map<int, std::vector<int>> &seedToTracksterAssociation) {
0152   // Protect from events with no seeding regions
0153   if (input.regions.empty())
0154     return;
0155 
0156   const int eventNumber = input.ev.eventAuxiliary().event();
0157   if (PatternRecognitionAlgoBaseT<TILES>::algo_verbosity_ > VerbosityLevel::Advanced) {
0158     edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "New Event";
0159   }
0160 
0161   edm::EventSetup const &es = input.es;
0162   const CaloGeometry &geom = es.getData(caloGeomToken_);
0163   rhtools_.setGeometry(geom);
0164 
0165   // Assume identical Z-positioning between positive and negative sides.
0166   // Also, layers inside the HGCAL geometry start from 1.
0167   for (unsigned int i = 0; i < rhtools_.lastLayer(); ++i) {
0168     layersPosZ_.push_back(rhtools_.getPositionLayer(i + 1).z());
0169     if (PatternRecognitionAlgoBaseT<TILES>::algo_verbosity_ > VerbosityLevel::Advanced) {
0170       edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "Layer " << i << " located at Z: " << layersPosZ_.back();
0171     }
0172   }
0173 
0174   clusters_.clear();
0175   tracksterSeedAlgoId_.clear();
0176 
0177   clusters_.resize(2 * rhtools_.lastLayer(false));
0178   std::vector<std::pair<int, int>> layerIdx2layerandSoa;  //used everywhere also to propagate cluster masking
0179 
0180   layerIdx2layerandSoa.reserve(input.layerClusters.size());
0181   unsigned int layerIdx = 0;
0182   for (auto const &lc : input.layerClusters) {
0183     if (input.mask[layerIdx] == 0.) {
0184       if (PatternRecognitionAlgoBaseT<TILES>::algo_verbosity_ > VerbosityLevel::Advanced) {
0185         edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "Skipping masked cluster: " << layerIdx;
0186       }
0187       layerIdx2layerandSoa.emplace_back(-1, -1);
0188       layerIdx++;
0189       continue;
0190     }
0191     const auto firstHitDetId = lc.hitsAndFractions()[0].first;
0192     int layer = rhtools_.getLayerWithOffset(firstHitDetId) - 1 +
0193                 rhtools_.lastLayer(false) * ((rhtools_.zside(firstHitDetId) + 1) >> 1);
0194     assert(layer >= 0);
0195     auto detId = lc.hitsAndFractions()[0].first;
0196     int layerClusterIndexInLayer = clusters_[layer].x.size();
0197     layerIdx2layerandSoa.emplace_back(layer, layerClusterIndexInLayer);
0198     float sum_x = 0.;
0199     float sum_y = 0.;
0200     float sum_sqr_x = 0.;
0201     float sum_sqr_y = 0.;
0202     float ref_x = lc.x();
0203     float ref_y = lc.y();
0204     float invClsize = 1. / lc.hitsAndFractions().size();
0205     for (auto const &hitsAndFractions : lc.hitsAndFractions()) {
0206       auto const &point = rhtools_.getPosition(hitsAndFractions.first);
0207       sum_x += point.x() - ref_x;
0208       sum_sqr_x += (point.x() - ref_x) * (point.x() - ref_x);
0209       sum_y += point.y() - ref_y;
0210       sum_sqr_y += (point.y() - ref_y) * (point.y() - ref_y);
0211     }
0212     // The variance of X for X uniform in circle of radius R is R^2/4,
0213     // therefore we multiply the sqrt(var) by 2 to have a rough estimate of the
0214     // radius. On the other hand, while averaging the x and y radius, we would
0215     // end up dividing by 2. Hence we omit the value here and in the average
0216     // below, too.
0217     float radius_x = sqrt((sum_sqr_x - (sum_x * sum_x) * invClsize) * invClsize);
0218     float radius_y = sqrt((sum_sqr_y - (sum_y * sum_y) * invClsize) * invClsize);
0219     if (PatternRecognitionAlgoBaseT<TILES>::algo_verbosity_ > VerbosityLevel::Advanced) {
0220       edm::LogVerbatim("PatternRecognitionbyCLUE3D")
0221           << "cluster rx: " << std::setw(5) << radius_x << ", ry: " << std::setw(5) << radius_y
0222           << ", r:  " << std::setw(5) << (radius_x + radius_y) << ", cells: " << std::setw(4)
0223           << lc.hitsAndFractions().size();
0224     }
0225 
0226     // The case of single cell layer clusters has to be handled differently.
0227 
0228     if (invClsize == 1.) {
0229       // Silicon case
0230       if (rhtools_.isSilicon(detId)) {
0231         radius_x = radius_y = rhtools_.getRadiusToSide(detId);
0232         if (PatternRecognitionAlgoBaseT<TILES>::algo_verbosity_ > VerbosityLevel::Advanced) {
0233           edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "Single cell cluster in silicon, rx: " << std::setw(5)
0234                                                          << radius_x << ", ry: " << std::setw(5) << radius_y;
0235         }
0236       } else {
0237         auto const &point = rhtools_.getPosition(detId);
0238         auto const &eta_phi_window = rhtools_.getScintDEtaDPhi(detId);
0239         radius_x = radius_y = point.perp() * eta_phi_window.second;
0240         if (PatternRecognitionAlgoBaseT<TILES>::algo_verbosity_ > VerbosityLevel::Advanced) {
0241           edm::LogVerbatim("PatternRecognitionbyCLUE3D")
0242               << "Single cell cluster in scintillator. rx: " << std::setw(5) << radius_x << ", ry: " << std::setw(5)
0243               << radius_y << ", eta-span: " << std::setw(5) << eta_phi_window.first << ", phi-span: " << std::setw(5)
0244               << eta_phi_window.second;
0245         }
0246       }
0247     }
0248 
0249     // Maybe check if these vectors can be reserved beforehands
0250     clusters_[layer].x.emplace_back(lc.x());
0251     clusters_[layer].y.emplace_back(lc.y());
0252     clusters_[layer].z.emplace_back(lc.z());
0253     clusters_[layer].r_over_absz.emplace_back(sqrt(lc.x() * lc.x() + lc.y() * lc.y()) / std::abs(lc.z()));
0254     clusters_[layer].radius.emplace_back(radius_x + radius_y);
0255     clusters_[layer].eta.emplace_back(lc.eta());
0256     clusters_[layer].phi.emplace_back(lc.phi());
0257     clusters_[layer].cells.push_back(lc.hitsAndFractions().size());
0258     clusters_[layer].algoId.push_back(lc.algo() - reco::CaloCluster::hgcal_em);
0259     clusters_[layer].isSilicon.push_back(rhtools_.isSilicon(detId));
0260     clusters_[layer].energy.emplace_back(lc.energy());
0261     clusters_[layer].isSeed.push_back(false);
0262     clusters_[layer].clusterIndex.emplace_back(-1);
0263     clusters_[layer].layerClusterOriginalIdx.emplace_back(layerIdx++);
0264     clusters_[layer].nearestHigher.emplace_back(-1, -1);
0265     clusters_[layer].rho.emplace_back(0.f);
0266     clusters_[layer].z_extension.emplace_back(0.f);
0267     clusters_[layer].delta.emplace_back(
0268         std::make_pair(std::numeric_limits<float>::max(), std::numeric_limits<int>::max()));
0269   }
0270   for (unsigned int layer = 0; layer < clusters_.size(); layer++) {
0271     clusters_[layer].followers.resize(clusters_[layer].x.size());
0272   }
0273 
0274   auto lastLayerPerSide = static_cast<int>(rhtools_.lastLayer(false));
0275   int maxLayer = 2 * lastLayerPerSide - 1;
0276   std::vector<int> numberOfClustersPerLayer(maxLayer, 0);
0277   for (int i = 0; i <= maxLayer; i++) {
0278     calculateLocalDensity(input.tiles, i, layerIdx2layerandSoa);
0279   }
0280   for (int i = 0; i <= maxLayer; i++) {
0281     calculateDistanceToHigher(input.tiles, i, layerIdx2layerandSoa);
0282   }
0283 
0284   auto nTracksters = findAndAssignTracksters(input.tiles, layerIdx2layerandSoa);
0285   if (PatternRecognitionAlgoBaseT<TILES>::algo_verbosity_ > VerbosityLevel::Advanced) {
0286     edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "Reconstructed " << nTracksters << " tracksters" << std::endl;
0287     dumpClusters(input.tiles, layerIdx2layerandSoa, eventNumber);
0288   }
0289 
0290   // Build Trackster
0291   result.resize(nTracksters);
0292 
0293   for (unsigned int layer = 0; layer < clusters_.size(); ++layer) {
0294     const auto &thisLayer = clusters_[layer];
0295     if (PatternRecognitionAlgoBaseT<TILES>::algo_verbosity_ > VerbosityLevel::Advanced) {
0296       edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "Examining Layer: " << layer;
0297     }
0298     for (unsigned int lc = 0; lc < thisLayer.x.size(); ++lc) {
0299       if (PatternRecognitionAlgoBaseT<TILES>::algo_verbosity_ > VerbosityLevel::Advanced) {
0300         edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "Trackster " << thisLayer.clusterIndex[lc];
0301       }
0302       if (thisLayer.clusterIndex[lc] >= 0) {
0303         if (PatternRecognitionAlgoBaseT<TILES>::algo_verbosity_ > VerbosityLevel::Advanced) {
0304           edm::LogVerbatim("PatternRecognitionbyCLUE3D") << " adding lcIdx: " << thisLayer.layerClusterOriginalIdx[lc];
0305         }
0306         result[thisLayer.clusterIndex[lc]].vertices().push_back(thisLayer.layerClusterOriginalIdx[lc]);
0307         result[thisLayer.clusterIndex[lc]].vertex_multiplicity().push_back(1);
0308         // loop over followers
0309         for (auto [follower_lyrIdx, follower_soaIdx] : thisLayer.followers[lc]) {
0310           std::array<unsigned int, 2> edge = {
0311               {(unsigned int)thisLayer.layerClusterOriginalIdx[lc],
0312                (unsigned int)clusters_[follower_lyrIdx].layerClusterOriginalIdx[follower_soaIdx]}};
0313           result[thisLayer.clusterIndex[lc]].edges().push_back(edge);
0314         }
0315       }
0316     }
0317   }
0318   size_t tracksterIndex = 0;
0319   result.erase(std::remove_if(std::begin(result),
0320                               std::end(result),
0321                               [&](auto const &v) {
0322                                 return static_cast<int>(v.vertices().size()) <
0323                                        minNumLayerCluster_[tracksterSeedAlgoId_[tracksterIndex++]];
0324                               }),
0325                result.end());
0326 
0327   result.shrink_to_fit();
0328 
0329   ticl::assignPCAtoTracksters(result,
0330                               input.layerClusters,
0331                               input.layerClustersTime,
0332                               rhtools_.getPositionLayer(rhtools_.lastLayerEE(false), false).z(),
0333                               rhtools_,
0334                               computeLocalTime_,
0335                               true,  // energy weighting
0336                               usePCACleaning_);
0337 
0338   if (PatternRecognitionAlgoBaseT<TILES>::algo_verbosity_ > VerbosityLevel::Advanced) {
0339     for (auto const &t : result) {
0340       edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "Barycenter: " << t.barycenter();
0341       edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "LCs: " << t.vertices().size();
0342       edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "Energy: " << t.raw_energy();
0343       edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "Regressed: " << t.regressed_energy();
0344     }
0345   }
0346 
0347   // Dump Tracksters information
0348   if (PatternRecognitionAlgoBaseT<TILES>::algo_verbosity_ > VerbosityLevel::Advanced) {
0349     dumpTracksters(layerIdx2layerandSoa, eventNumber, result);
0350   }
0351 
0352   // Reset internal clusters_ structure of array for next event
0353   reset();
0354   if (PatternRecognitionAlgoBaseT<TILES>::algo_verbosity_ > VerbosityLevel::Advanced) {
0355     edm::LogVerbatim("PatternRecognitionbyCLUE3D") << std::endl;
0356   }
0357 }
0358 template <typename TILES>
0359 void PatternRecognitionbyCLUE3D<TILES>::filter(std::vector<Trackster> &output,
0360                                                const std::vector<Trackster> &inTracksters,
0361                                                const typename PatternRecognitionAlgoBaseT<TILES>::Inputs &input,
0362                                                std::unordered_map<int, std::vector<int>> &seedToTracksterAssociation) {
0363   auto isHAD = [this](const Trackster &t) -> bool {
0364     auto const hadProb = t.id_probability(ticl::Trackster::ParticleType::charged_hadron) +
0365                          t.id_probability(ticl::Trackster::ParticleType::neutral_hadron);
0366     return hadProb >= cutHadProb_;
0367   };
0368 
0369   if (doPidCut_) {
0370     for (auto const &t : inTracksters) {
0371       if (!isHAD(t)) {
0372         output.push_back(t);
0373       }
0374     }
0375   } else {
0376     output = inTracksters;
0377   }
0378 }
0379 
0380 template <typename TILES>
0381 void PatternRecognitionbyCLUE3D<TILES>::calculateLocalDensity(
0382     const TILES &tiles, const int layerId, const std::vector<std::pair<int, int>> &layerIdx2layerandSoa) {
0383   constexpr int nEtaBin = TILES::constants_type_t::nEtaBins;
0384   constexpr int nPhiBin = TILES::constants_type_t::nPhiBins;
0385   auto &clustersOnLayer = clusters_[layerId];
0386   unsigned int numberOfClusters = clustersOnLayer.x.size();
0387 
0388   auto isReachable = [](float r0, float r1, float phi0, float phi1, float delta_sqr) -> bool {
0389     auto delta_phi = reco::deltaPhi(phi0, phi1);
0390     return (r0 - r1) * (r0 - r1) + r1 * r1 * delta_phi * delta_phi < delta_sqr;
0391   };
0392   auto distance_debug = [&](float x1, float x2, float y1, float y2) -> float {
0393     return sqrt((x1 - x2) * (x1 - x2) + (y1 - y2) * (y1 - y2));
0394   };
0395 
0396   for (unsigned int i = 0; i < numberOfClusters; i++) {
0397     auto algoId = clustersOnLayer.algoId[i];
0398     // We need to partition the two sides of the HGCAL detector
0399     auto lastLayerPerSide = static_cast<int>(rhtools_.lastLayer(false));
0400     int minLayer = 0;
0401     int maxLayer = 2 * lastLayerPerSide - 1;
0402     if (layerId < lastLayerPerSide) {
0403       minLayer = std::max(layerId - densitySiblingLayers_[algoId], minLayer);
0404       maxLayer = std::min(layerId + densitySiblingLayers_[algoId], lastLayerPerSide - 1);
0405     } else {
0406       minLayer = std::max(layerId - densitySiblingLayers_[algoId], lastLayerPerSide);
0407       maxLayer = std::min(layerId + densitySiblingLayers_[algoId], maxLayer);
0408     }
0409     float deltaLayersZ = std::abs(layersPosZ_[maxLayer % lastLayerPerSide] - layersPosZ_[minLayer % lastLayerPerSide]);
0410 
0411     for (int currentLayer = minLayer; currentLayer <= maxLayer; currentLayer++) {
0412       if (PatternRecognitionAlgoBaseT<TILES>::algo_verbosity_ > VerbosityLevel::Advanced) {
0413         edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "RefLayer: " << layerId << " SoaIDX: " << i;
0414         edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "NextLayer: " << currentLayer;
0415       }
0416       const auto &tileOnLayer = tiles[currentLayer];
0417       bool onSameLayer = (currentLayer == layerId);
0418       if (PatternRecognitionAlgoBaseT<TILES>::algo_verbosity_ > VerbosityLevel::Advanced) {
0419         edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "onSameLayer: " << onSameLayer;
0420       }
0421       const int etaWindow = 2;
0422       const int phiWindow = 2;
0423       int etaBinMin = std::max(tileOnLayer.etaBin(clustersOnLayer.eta[i]) - etaWindow, 0);
0424       int etaBinMax = std::min(tileOnLayer.etaBin(clustersOnLayer.eta[i]) + etaWindow, nEtaBin - 1);
0425       int phiBinMin = tileOnLayer.phiBin(clustersOnLayer.phi[i]) - phiWindow;
0426       int phiBinMax = tileOnLayer.phiBin(clustersOnLayer.phi[i]) + phiWindow;
0427       if (PatternRecognitionAlgoBaseT<TILES>::algo_verbosity_ > VerbosityLevel::Advanced) {
0428         edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "eta: " << clustersOnLayer.eta[i];
0429         edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "phi: " << clustersOnLayer.phi[i];
0430         edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "etaBinMin: " << etaBinMin << ", etaBinMax: " << etaBinMax;
0431         edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "phiBinMin: " << phiBinMin << ", phiBinMax: " << phiBinMax;
0432       }
0433       for (int ieta = etaBinMin; ieta <= etaBinMax; ++ieta) {
0434         auto offset = ieta * nPhiBin;
0435         if (PatternRecognitionAlgoBaseT<TILES>::algo_verbosity_ > VerbosityLevel::Advanced) {
0436           edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "offset: " << offset;
0437         }
0438         for (int iphi_it = phiBinMin; iphi_it <= phiBinMax; ++iphi_it) {
0439           int iphi = ((iphi_it % nPhiBin + nPhiBin) % nPhiBin);
0440           if (PatternRecognitionAlgoBaseT<TILES>::algo_verbosity_ > VerbosityLevel::Advanced) {
0441             edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "iphi: " << iphi;
0442             edm::LogVerbatim("PatternRecognitionbyCLUE3D")
0443                 << "Entries in tileBin: " << tileOnLayer[offset + iphi].size();
0444           }
0445           for (auto otherClusterIdx : tileOnLayer[offset + iphi]) {
0446             auto const &layerandSoa = layerIdx2layerandSoa[otherClusterIdx];
0447             // Skip masked layer clusters
0448             if ((layerandSoa.first == -1) && (layerandSoa.second == -1)) {
0449               if (PatternRecognitionAlgoBaseT<TILES>::algo_verbosity_ > VerbosityLevel::Advanced) {
0450                 edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "Skipping masked layerIdx " << otherClusterIdx;
0451               }
0452               continue;
0453             }
0454             auto const &clustersLayer = clusters_[layerandSoa.first];
0455             if (PatternRecognitionAlgoBaseT<TILES>::algo_verbosity_ > VerbosityLevel::Advanced) {
0456               edm::LogVerbatim("PatternRecognitionbyCLUE3D")
0457                   << "OtherLayer: " << layerandSoa.first << " SoaIDX: " << layerandSoa.second;
0458               edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "OtherEta: " << clustersLayer.eta[layerandSoa.second];
0459               edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "OtherPhi: " << clustersLayer.phi[layerandSoa.second];
0460             }
0461 
0462             bool onSameCluster = clustersOnLayer.layerClusterOriginalIdx[i] == otherClusterIdx;
0463             if (onSameLayer && !densityOnSameLayer_ && !onSameCluster) {
0464               if (PatternRecognitionAlgoBaseT<TILES>::algo_verbosity_ > VerbosityLevel::Advanced) {
0465                 edm::LogVerbatim("PatternRecognitionbyCLUE3D")
0466                     << "Skipping different cluster " << otherClusterIdx << "in the same layer " << currentLayer;
0467               }
0468               continue;
0469             }
0470 
0471             bool reachable = false;
0472             if (useAbsoluteProjectiveScale_) {
0473               if (useClusterDimensionXY_) {
0474                 reachable = isReachable(clustersOnLayer.r_over_absz[i] * clustersOnLayer.z[i],
0475                                         clustersLayer.r_over_absz[layerandSoa.second] * clustersOnLayer.z[i],
0476                                         clustersOnLayer.phi[i],
0477                                         clustersLayer.phi[layerandSoa.second],
0478                                         clustersOnLayer.radius[i] * clustersOnLayer.radius[i]);
0479               } else {
0480                 // Still differentiate between silicon and Scintillator.
0481                 // Scintillator has yet to be studied further.
0482                 if (clustersOnLayer.isSilicon[i]) {
0483                   reachable = isReachable(clustersOnLayer.r_over_absz[i] * clustersOnLayer.z[i],
0484                                           clustersLayer.r_over_absz[layerandSoa.second] * clustersOnLayer.z[i],
0485                                           clustersOnLayer.phi[i],
0486                                           clustersLayer.phi[layerandSoa.second],
0487                                           densityXYDistanceSqr_[algoId]);
0488                 } else {
0489                   reachable = isReachable(clustersOnLayer.r_over_absz[i] * clustersOnLayer.z[i],
0490                                           clustersLayer.r_over_absz[layerandSoa.second] * clustersOnLayer.z[i],
0491                                           clustersOnLayer.phi[i],
0492                                           clustersLayer.phi[layerandSoa.second],
0493                                           clustersOnLayer.radius[i] * clustersOnLayer.radius[i]);
0494                 }
0495               }
0496             } else {
0497               reachable = (reco::deltaR2(clustersOnLayer.eta[i],
0498                                          clustersOnLayer.phi[i],
0499                                          clustersLayer.eta[layerandSoa.second],
0500                                          clustersLayer.phi[layerandSoa.second]) < densityEtaPhiDistanceSqr_[algoId]);
0501             }
0502             if (PatternRecognitionAlgoBaseT<TILES>::algo_verbosity_ > VerbosityLevel::Advanced) {
0503               edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "Distance[eta,phi]: "
0504                                                              << reco::deltaR2(clustersOnLayer.eta[i],
0505                                                                               clustersOnLayer.phi[i],
0506                                                                               clustersLayer.eta[layerandSoa.second],
0507                                                                               clustersLayer.phi[layerandSoa.second]);
0508               auto dist = distance_debug(
0509                   clustersOnLayer.r_over_absz[i],
0510                   clustersLayer.r_over_absz[layerandSoa.second],
0511                   clustersOnLayer.r_over_absz[i] * std::abs(clustersOnLayer.phi[i]),
0512                   clustersLayer.r_over_absz[layerandSoa.second] * std::abs(clustersLayer.phi[layerandSoa.second]));
0513               edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "Distance[cm]: " << (dist * clustersOnLayer.z[i]);
0514               edm::LogVerbatim("PatternRecognitionbyCLUE3D")
0515                   << "Energy Other:   " << clustersLayer.energy[layerandSoa.second];
0516               edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "Cluster radius: " << clustersOnLayer.radius[i];
0517             }
0518             if (reachable) {
0519               float factor_same_layer_different_cluster = (onSameLayer && !densityOnSameLayer_) ? 0.f : 1.f;
0520               auto energyToAdd = (clustersOnLayer.layerClusterOriginalIdx[i] == otherClusterIdx
0521                                       ? 1.f
0522                                       : kernelDensityFactor_[algoId] * factor_same_layer_different_cluster) *
0523                                  clustersLayer.energy[layerandSoa.second];
0524               clustersOnLayer.rho[i] += energyToAdd;
0525               clustersOnLayer.z_extension[i] = deltaLayersZ;
0526               if (PatternRecognitionAlgoBaseT<TILES>::algo_verbosity_ > VerbosityLevel::Advanced) {
0527                 edm::LogVerbatim("PatternRecognitionbyCLUE3D")
0528                     << "Adding " << energyToAdd << " partial " << clustersOnLayer.rho[i];
0529               }
0530             }
0531           }  // end of loop on possible compatible clusters
0532         }  // end of loop over phi-bin region
0533       }  // end of loop over eta-bin region
0534     }  // end of loop on the sibling layers
0535     if (rescaleDensityByZ_) {
0536       if (PatternRecognitionAlgoBaseT<TILES>::algo_verbosity_ > VerbosityLevel::Advanced) {
0537         edm::LogVerbatim("PatternRecognitionbyCLUE3D")
0538             << "Rescaling original density: " << clustersOnLayer.rho[i] << " by Z: " << deltaLayersZ
0539             << " to final density/cm: " << clustersOnLayer.rho[i] / deltaLayersZ;
0540       }
0541       clustersOnLayer.rho[i] /= deltaLayersZ;
0542     }
0543   }  // end of loop over clusters on this layer
0544 }
0545 
0546 template <typename TILES>
0547 void PatternRecognitionbyCLUE3D<TILES>::calculateDistanceToHigher(
0548     const TILES &tiles, const int layerId, const std::vector<std::pair<int, int>> &layerIdx2layerandSoa) {
0549   constexpr int nEtaBin = TILES::constants_type_t::nEtaBins;
0550   constexpr int nPhiBin = TILES::constants_type_t::nPhiBins;
0551   auto &clustersOnLayer = clusters_[layerId];
0552   unsigned int numberOfClusters = clustersOnLayer.x.size();
0553 
0554   auto distanceSqr = [](float r0, float r1, float phi0, float phi1) -> float {
0555     auto delta_phi = reco::deltaPhi(phi0, phi1);
0556     return (r0 - r1) * (r0 - r1) + r1 * r1 * delta_phi * delta_phi;
0557   };
0558 
0559   for (unsigned int i = 0; i < numberOfClusters; i++) {
0560     if (PatternRecognitionAlgoBaseT<TILES>::algo_verbosity_ > VerbosityLevel::Advanced) {
0561       edm::LogVerbatim("PatternRecognitionbyCLUE3D")
0562           << "Starting searching nearestHigher on " << layerId << " with rho: " << clustersOnLayer.rho[i]
0563           << " at eta, phi: " << tiles[layerId].etaBin(clustersOnLayer.eta[i]) << ", "
0564           << tiles[layerId].phiBin(clustersOnLayer.phi[i]);
0565     }
0566     // We need to partition the two sides of the HGCAL detector
0567     auto lastLayerPerSide = static_cast<int>(rhtools_.lastLayer(false));
0568     int minLayer = 0;
0569     int maxLayer = 2 * lastLayerPerSide - 1;
0570     auto algoId = clustersOnLayer.algoId[i];
0571     if (layerId < lastLayerPerSide) {
0572       minLayer = std::max(layerId - densitySiblingLayers_[algoId], minLayer);
0573       maxLayer = std::min(layerId + densitySiblingLayers_[algoId], lastLayerPerSide - 1);
0574     } else {
0575       minLayer = std::max(layerId - densitySiblingLayers_[algoId], lastLayerPerSide + 1);
0576       maxLayer = std::min(layerId + densitySiblingLayers_[algoId], maxLayer);
0577     }
0578     constexpr float maxDelta = std::numeric_limits<float>::max();
0579     float i_delta = maxDelta;
0580     std::pair<int, int> i_nearestHigher(-1, -1);
0581     std::pair<float, int> nearest_distances(maxDelta, std::numeric_limits<int>::max());
0582     for (int currentLayer = minLayer; currentLayer <= maxLayer; currentLayer++) {
0583       if (!nearestHigherOnSameLayer_ && (layerId == currentLayer))
0584         continue;
0585       const auto &tileOnLayer = tiles[currentLayer];
0586       int etaWindow = 1;
0587       int phiWindow = 1;
0588       int etaBinMin = std::max(tileOnLayer.etaBin(clustersOnLayer.eta[i]) - etaWindow, 0);
0589       int etaBinMax = std::min(tileOnLayer.etaBin(clustersOnLayer.eta[i]) + etaWindow, nEtaBin - 1);
0590       int phiBinMin = tileOnLayer.phiBin(clustersOnLayer.phi[i]) - phiWindow;
0591       int phiBinMax = tileOnLayer.phiBin(clustersOnLayer.phi[i]) + phiWindow;
0592       for (int ieta = etaBinMin; ieta <= etaBinMax; ++ieta) {
0593         auto offset = ieta * nPhiBin;
0594         for (int iphi_it = phiBinMin; iphi_it <= phiBinMax; ++iphi_it) {
0595           int iphi = ((iphi_it % nPhiBin + nPhiBin) % nPhiBin);
0596           if (PatternRecognitionAlgoBaseT<TILES>::algo_verbosity_ > VerbosityLevel::Advanced) {
0597             edm::LogVerbatim("PatternRecognitionbyCLUE3D")
0598                 << "Searching nearestHigher on " << currentLayer << " eta, phi: " << ieta << ", " << iphi_it << " "
0599                 << iphi << " " << offset << " " << (offset + iphi);
0600           }
0601           for (auto otherClusterIdx : tileOnLayer[offset + iphi]) {
0602             auto const &layerandSoa = layerIdx2layerandSoa[otherClusterIdx];
0603             // Skip masked layer clusters
0604             if ((layerandSoa.first == -1) && (layerandSoa.second == -1))
0605               continue;
0606             auto const &clustersOnOtherLayer = clusters_[layerandSoa.first];
0607             auto dist = maxDelta;
0608             auto dist_transverse = maxDelta;
0609             int dist_layers = std::abs(layerandSoa.first - layerId);
0610             if (useAbsoluteProjectiveScale_) {
0611               dist_transverse = distanceSqr(clustersOnLayer.r_over_absz[i] * clustersOnLayer.z[i],
0612                                             clustersOnOtherLayer.r_over_absz[layerandSoa.second] * clustersOnLayer.z[i],
0613                                             clustersOnLayer.phi[i],
0614                                             clustersOnOtherLayer.phi[layerandSoa.second]);
0615               // Add Z-scale to the final distance
0616               dist = dist_transverse;
0617             } else {
0618               dist = reco::deltaR2(clustersOnLayer.eta[i],
0619                                    clustersOnLayer.phi[i],
0620                                    clustersOnOtherLayer.eta[layerandSoa.second],
0621                                    clustersOnOtherLayer.phi[layerandSoa.second]);
0622               dist_transverse = dist;
0623             }
0624             bool foundHigher = (clustersOnOtherLayer.rho[layerandSoa.second] > clustersOnLayer.rho[i]) ||
0625                                (clustersOnOtherLayer.rho[layerandSoa.second] == clustersOnLayer.rho[i] &&
0626                                 clustersOnOtherLayer.layerClusterOriginalIdx[layerandSoa.second] >
0627                                     clustersOnLayer.layerClusterOriginalIdx[i]);
0628             if (PatternRecognitionAlgoBaseT<TILES>::algo_verbosity_ > VerbosityLevel::Advanced) {
0629               edm::LogVerbatim("PatternRecognitionbyCLUE3D")
0630                   << "Searching nearestHigher on " << currentLayer
0631                   << " with rho: " << clustersOnOtherLayer.rho[layerandSoa.second]
0632                   << " on layerIdxInSOA: " << layerandSoa.first << ", " << layerandSoa.second
0633                   << " with distance: " << sqrt(dist) << " foundHigher: " << foundHigher;
0634             }
0635             if (foundHigher && dist <= i_delta) {
0636               // update i_delta
0637               i_delta = dist;
0638               nearest_distances = std::make_pair(sqrt(dist_transverse), dist_layers);
0639               // update i_nearestHigher
0640               i_nearestHigher = layerandSoa;
0641             }
0642           }  // End of loop on clusters
0643         }  // End of loop on phi bins
0644       }  // End of loop on eta bins
0645     }  // End of loop on layers
0646 
0647     clustersOnLayer.delta[i] = nearest_distances;
0648     clustersOnLayer.nearestHigher[i] = i_nearestHigher;
0649   }  // End of loop on clusters
0650 }
0651 
0652 template <typename TILES>
0653 int PatternRecognitionbyCLUE3D<TILES>::findAndAssignTracksters(
0654     const TILES &tiles, const std::vector<std::pair<int, int>> &layerIdx2layerandSoa) {
0655   unsigned int nTracksters = 0;
0656   std::vector<std::pair<int, int>> localStack;
0657   const auto &critical_transverse_distance =
0658       useAbsoluteProjectiveScale_ ? criticalXYDistance_ : criticalEtaPhiDistance_;
0659   // find cluster seeds and outlier
0660   for (unsigned int layer = 0; layer < 2 * rhtools_.lastLayer(); layer++) {
0661     auto &clustersOnLayer = clusters_[layer];
0662     unsigned int numberOfClusters = clustersOnLayer.x.size();
0663     for (unsigned int i = 0; i < numberOfClusters; i++) {
0664       // initialize clusterIndex
0665       clustersOnLayer.clusterIndex[i] = -1;
0666       auto algoId = clustersOnLayer.algoId[i];
0667       bool isSeed = (clustersOnLayer.delta[i].first > critical_transverse_distance[algoId] ||
0668                      clustersOnLayer.delta[i].second > criticalZDistanceLyr_[algoId]) &&
0669                     (clustersOnLayer.rho[i] >= criticalDensity_[algoId]) &&
0670                     (clustersOnLayer.energy[i] / clustersOnLayer.rho[i] > criticalSelfDensity_[algoId]);
0671       if (!clustersOnLayer.isSilicon[i]) {
0672         isSeed = (clustersOnLayer.delta[i].first > clustersOnLayer.radius[i] ||
0673                   clustersOnLayer.delta[i].second > criticalZDistanceLyr_[algoId]) &&
0674                  (clustersOnLayer.rho[i] >= criticalDensity_[algoId]) &&
0675                  (clustersOnLayer.energy[i] / clustersOnLayer.rho[i] > criticalSelfDensity_[algoId]);
0676       }
0677       bool isOutlier =
0678           (clustersOnLayer.delta[i].first > outlierMultiplier_[algoId] * critical_transverse_distance[algoId]) &&
0679           (clustersOnLayer.rho[i] < criticalDensity_[algoId]);
0680       if (isSeed) {
0681         if (PatternRecognitionAlgoBaseT<TILES>::algo_verbosity_ > VerbosityLevel::Advanced) {
0682           edm::LogVerbatim("PatternRecognitionbyCLUE3D")
0683               << "Found seed on Layer " << layer << " SOAidx: " << i << " assigned ClusterIdx: " << nTracksters;
0684         }
0685         clustersOnLayer.clusterIndex[i] = nTracksters++;
0686         tracksterSeedAlgoId_.push_back(algoId);
0687         clustersOnLayer.isSeed[i] = true;
0688         localStack.emplace_back(layer, i);
0689       } else if (!isOutlier) {
0690         auto [lyrIdx, soaIdx] = clustersOnLayer.nearestHigher[i];
0691         if (PatternRecognitionAlgoBaseT<TILES>::algo_verbosity_ > VerbosityLevel::Advanced) {
0692           edm::LogVerbatim("PatternRecognitionbyCLUE3D")
0693               << "Found follower on Layer " << layer << " SOAidx: " << i << " attached to cluster on layer: " << lyrIdx
0694               << " SOAidx: " << soaIdx;
0695         }
0696         if (lyrIdx >= 0)
0697           clusters_[lyrIdx].followers[soaIdx].emplace_back(layer, i);
0698       } else {
0699         if (PatternRecognitionAlgoBaseT<TILES>::algo_verbosity_ > VerbosityLevel::Advanced) {
0700           edm::LogVerbatim("PatternRecognitionbyCLUE3D")
0701               << "Found Outlier on Layer " << layer << " SOAidx: " << i << " with rho: " << clustersOnLayer.rho[i]
0702               << " and delta: " << clustersOnLayer.delta[i].first << ", " << clustersOnLayer.delta[i].second;
0703         }
0704       }
0705     }
0706   }
0707 
0708   // Propagate cluster index
0709   while (!localStack.empty()) {
0710     auto [lyrIdx, soaIdx] = localStack.back();
0711     auto &thisSeed = clusters_[lyrIdx].followers[soaIdx];
0712     localStack.pop_back();
0713 
0714     // loop over followers
0715     for (auto [follower_lyrIdx, follower_soaIdx] : thisSeed) {
0716       // pass id to a follower
0717       clusters_[follower_lyrIdx].clusterIndex[follower_soaIdx] = clusters_[lyrIdx].clusterIndex[soaIdx];
0718       // push this follower to localStack
0719       localStack.emplace_back(follower_lyrIdx, follower_soaIdx);
0720     }
0721   }
0722   return nTracksters;
0723 }
0724 
0725 template <typename TILES>
0726 void PatternRecognitionbyCLUE3D<TILES>::fillPSetDescription(edm::ParameterSetDescription &iDesc) {
0727   iDesc.add<int>("algo_verbosity", 0);
0728   iDesc.add<std::vector<double>>("criticalDensity", {4, 4, 4})->setComment("in GeV");
0729   iDesc.add<std::vector<double>>("criticalSelfDensity", {0.15, 0.15, 0.15} /* roughly 1/(densitySiblingLayers+1) */)
0730       ->setComment("Minimum ratio of self_energy/local_density to become a seed.");
0731   iDesc.add<std::vector<int>>("densitySiblingLayers", {3, 3, 3})
0732       ->setComment(
0733           "inclusive, layers to consider while computing local density and searching for nearestHigher higher");
0734   iDesc.add<std::vector<double>>("densityEtaPhiDistanceSqr", {0.0008, 0.0008, 0.0008})
0735       ->setComment("in eta,phi space, distance to consider for local density");
0736   iDesc.add<std::vector<double>>("densityXYDistanceSqr", {3.24, 3.24, 3.24})
0737       ->setComment("in cm, distance on the transverse plane to consider for local density");
0738   iDesc.add<std::vector<double>>("kernelDensityFactor", {0.2, 0.2, 0.2})
0739       ->setComment("Kernel factor to be applied to other LC while computing the local density");
0740   iDesc.add<bool>("densityOnSameLayer", false);
0741   iDesc.add<bool>("nearestHigherOnSameLayer", false)
0742       ->setComment("Allow the nearestHigher to be located on the same layer");
0743   iDesc.add<bool>("useAbsoluteProjectiveScale", true)
0744       ->setComment("Express all cuts in terms of r/z*z_0{,phi} projective variables");
0745   iDesc.add<bool>("useClusterDimensionXY", false)
0746       ->setComment(
0747           "Boolean. If true use the estimated cluster radius to determine the cluster compatibility while computing "
0748           "the local density");
0749   iDesc.add<bool>("rescaleDensityByZ", false)
0750       ->setComment(
0751           "Rescale local density by the extension of the Z 'volume' explored. The transvere dimension is, at present, "
0752           "fixed and factored out.");
0753   iDesc.add<std::vector<double>>("criticalEtaPhiDistance", {0.025, 0.025, 0.025})
0754       ->setComment("Minimal distance in eta,phi space from nearestHigher to become a seed");
0755   iDesc.add<std::vector<double>>("criticalXYDistance", {1.8, 1.8, 1.8})
0756       ->setComment("Minimal distance in cm on the XY plane from nearestHigher to become a seed");
0757   iDesc.add<std::vector<int>>("criticalZDistanceLyr", {5, 5, 5})
0758       ->setComment("Minimal distance in layers along the Z axis from nearestHigher to become a seed");
0759   iDesc.add<std::vector<double>>("outlierMultiplier", {2, 2, 2})
0760       ->setComment("Minimal distance in transverse space from nearestHigher to become an outlier");
0761   iDesc.add<std::vector<int>>("minNumLayerCluster", {2, 2, 2})->setComment("Not Inclusive");
0762   iDesc.add<bool>("doPidCut", false);
0763   iDesc.add<double>("cutHadProb", 0.5);
0764   iDesc.add<bool>("computeLocalTime", false);
0765   iDesc.add<bool>("usePCACleaning", false)->setComment("Enable PCA cleaning alorithm");
0766 }
0767 
0768 template class ticl::PatternRecognitionbyCLUE3D<TICLLayerTiles>;
0769 template class ticl::PatternRecognitionbyCLUE3D<TICLLayerTilesHFNose>;