Back to home page

Project CMSSW displayed by LXR
 
 

    

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

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       eidInputName_(conf.getParameter<std::string>("eid_input_name")),
0043       eidOutputNameEnergy_(conf.getParameter<std::string>("eid_output_name_energy")),
0044       eidOutputNameId_(conf.getParameter<std::string>("eid_output_name_id")),
0045       eidMinClusterEnergy_(conf.getParameter<double>("eid_min_cluster_energy")),
0046       eidNLayers_(conf.getParameter<int>("eid_n_layers")),
0047       eidNClusters_(conf.getParameter<int>("eid_n_clusters")),
0048       computeLocalTime_(conf.getParameter<bool>("computeLocalTime")){};
0049 
0050 template <typename TILES>
0051 void PatternRecognitionbyCLUE3D<TILES>::dumpTiles(const TILES &tiles) const {
0052   constexpr int nEtaBin = TILES::constants_type_t::nEtaBins;
0053   constexpr int nPhiBin = TILES::constants_type_t::nPhiBins;
0054   auto lastLayerPerSide = static_cast<int>(rhtools_.lastLayer(false));
0055   int maxLayer = 2 * lastLayerPerSide - 1;
0056   for (int layer = 0; layer <= maxLayer; layer++) {
0057     for (int ieta = 0; ieta < nEtaBin; ieta++) {
0058       auto offset = ieta * nPhiBin;
0059       for (int phi = 0; phi < nPhiBin; phi++) {
0060         int iphi = ((phi % nPhiBin + nPhiBin) % nPhiBin);
0061         if (!tiles[layer][offset + iphi].empty()) {
0062           if (this->algo_verbosity_ > VerbosityLevel::Advanced) {
0063             edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "Layer: " << layer << " ieta: " << ieta << " phi: " << phi
0064                                                            << " " << tiles[layer][offset + iphi].size();
0065           }
0066         }
0067       }
0068     }
0069   }
0070 }
0071 
0072 template <typename TILES>
0073 void PatternRecognitionbyCLUE3D<TILES>::dumpTracksters(const std::vector<std::pair<int, int>> &layerIdx2layerandSoa,
0074                                                        const int eventNumber,
0075                                                        const std::vector<Trackster> &tracksters) const {
0076   if (PatternRecognitionAlgoBaseT<TILES>::algo_verbosity_ > VerbosityLevel::Advanced) {
0077     edm::LogVerbatim("PatternRecognitionbyCLUE3D")
0078         << "[evt, tracksterId, cells, prob_photon, prob_ele, prob_chad, prob_nhad, layer_i, x_i, y_i, eta_i, phi_i, "
0079            "energy_i, radius_i, rho_i, z_extension, delta_tr, delta_lyr, isSeed_i";
0080   }
0081 
0082   int num = 0;
0083   const std::string sep(", ");
0084   for (auto const &t : tracksters) {
0085     for (auto v : t.vertices()) {
0086       auto [lyrIdx, soaIdx] = layerIdx2layerandSoa[v];
0087       auto const &thisLayer = clusters_[lyrIdx];
0088       if (PatternRecognitionAlgoBaseT<TILES>::algo_verbosity_ > VerbosityLevel::Advanced) {
0089         edm::LogVerbatim("PatternRecognitionbyCLUE3D_NTP")
0090             << std::setw(4) << eventNumber << sep << std::setw(4) << num << sep << std::setw(4) << t.vertices().size()
0091             << sep << std::setw(8) << t.id_probability(ticl::Trackster::ParticleType::photon) << sep << std::setw(8)
0092             << t.id_probability(ticl::Trackster::ParticleType::electron) << sep << std::setw(8)
0093             << t.id_probability(ticl::Trackster::ParticleType::charged_hadron) << sep << std::setw(8)
0094             << t.id_probability(ticl::Trackster::ParticleType::neutral_hadron) << sep << std::setw(4) << lyrIdx << sep
0095             << std::setw(10) << thisLayer.x[soaIdx] << sep << std::setw(10) << thisLayer.y[soaIdx] << sep
0096             << std::setw(10) << thisLayer.eta[soaIdx] << sep << std::setw(10) << thisLayer.phi[soaIdx] << sep
0097             << std::setw(10) << thisLayer.energy[soaIdx] << sep << std::setw(10) << thisLayer.radius[soaIdx] << sep
0098             << std::setw(10) << thisLayer.rho[soaIdx] << sep << std::setw(10) << thisLayer.z_extension[soaIdx] << sep
0099             << std::setw(10) << thisLayer.delta[soaIdx].first << sep << std::setw(10) << thisLayer.delta[soaIdx].second
0100             << sep << std::setw(4) << thisLayer.isSeed[soaIdx];
0101       }
0102     }
0103     num++;
0104   }
0105 }
0106 
0107 template <typename TILES>
0108 void PatternRecognitionbyCLUE3D<TILES>::dumpClusters(const TILES &tiles,
0109                                                      const std::vector<std::pair<int, int>> &layerIdx2layerandSoa,
0110                                                      const int eventNumber) const {
0111   if (PatternRecognitionAlgoBaseT<TILES>::algo_verbosity_ > VerbosityLevel::Advanced) {
0112     edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "[evt, lyr, Seed,      x,       y,       z, r/|z|,   eta,   phi, "
0113                                                       "etab,  phib, cells, enrgy, e/rho,   rho,   z_ext, "
0114                                                       "   dlt_tr,   dlt_lyr, "
0115                                                       " nestHL, nestHSoaIdx, radius, clIdx, lClOrigIdx, SOAidx";
0116   }
0117 
0118   for (unsigned int layer = 0; layer < clusters_.size(); layer++) {
0119     auto const &thisLayer = clusters_[layer];
0120     int num = 0;
0121     for (auto v : thisLayer.x) {
0122       if (PatternRecognitionAlgoBaseT<TILES>::algo_verbosity_ > VerbosityLevel::Advanced) {
0123         edm::LogVerbatim("PatternRecognitionbyCLUE3D")
0124             << std::setw(4) << eventNumber << ", " << std::setw(3) << layer << ", " << std::setw(4)
0125             << thisLayer.isSeed[num] << ", " << std::setprecision(3) << std::fixed << v << ", " << thisLayer.y[num]
0126             << ", " << thisLayer.z[num] << ", " << thisLayer.r_over_absz[num] << ", " << thisLayer.eta[num] << ", "
0127             << thisLayer.phi[num] << ", " << std::setw(5) << tiles[layer].etaBin(thisLayer.eta[num]) << ", "
0128             << std::setw(5) << tiles[layer].phiBin(thisLayer.phi[num]) << ", " << std::setw(4) << thisLayer.cells[num]
0129             << ", " << std::setprecision(3) << thisLayer.energy[num] << ", "
0130             << (thisLayer.energy[num] / thisLayer.rho[num]) << ", " << thisLayer.rho[num] << ", "
0131             << thisLayer.z_extension[num] << ", " << std::scientific << thisLayer.delta[num].first << ", "
0132             << std::setw(10) << thisLayer.delta[num].second << ", " << std::setw(5)
0133             << thisLayer.nearestHigher[num].first << ", " << std::setw(10) << thisLayer.nearestHigher[num].second
0134             << ", " << std::defaultfloat << std::setprecision(3) << thisLayer.radius[num] << ", " << std::setw(5)
0135             << thisLayer.clusterIndex[num] << ", " << std::setw(4) << thisLayer.layerClusterOriginalIdx[num] << ", "
0136             << std::setw(4) << num << ", ClusterInfo";
0137       }
0138       ++num;
0139     }
0140   }
0141   for (unsigned int lcIdx = 0; lcIdx < layerIdx2layerandSoa.size(); lcIdx++) {
0142     auto const &layerandSoa = layerIdx2layerandSoa[lcIdx];
0143     // Skip masked layer clusters
0144     if ((layerandSoa.first == -1) && (layerandSoa.second == -1))
0145       continue;
0146     if (PatternRecognitionAlgoBaseT<TILES>::algo_verbosity_ > VerbosityLevel::Advanced) {
0147       edm::LogVerbatim("PatternRecognitionbyCLUE3D")
0148           << "lcIdx: " << lcIdx << " on Layer: " << layerandSoa.first << " SOA: " << layerandSoa.second;
0149     }
0150   }
0151 }
0152 
0153 template <typename TILES>
0154 void PatternRecognitionbyCLUE3D<TILES>::makeTracksters(
0155     const typename PatternRecognitionAlgoBaseT<TILES>::Inputs &input,
0156     std::vector<Trackster> &result,
0157     std::unordered_map<int, std::vector<int>> &seedToTracksterAssociation) {
0158   // Protect from events with no seeding regions
0159   if (input.regions.empty())
0160     return;
0161 
0162   const int eventNumber = input.ev.eventAuxiliary().event();
0163   if (PatternRecognitionAlgoBaseT<TILES>::algo_verbosity_ > VerbosityLevel::Advanced) {
0164     edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "New Event";
0165   }
0166 
0167   edm::EventSetup const &es = input.es;
0168   const CaloGeometry &geom = es.getData(caloGeomToken_);
0169   rhtools_.setGeometry(geom);
0170 
0171   // Assume identical Z-positioning between positive and negative sides.
0172   // Also, layers inside the HGCAL geometry start from 1.
0173   for (unsigned int i = 0; i < rhtools_.lastLayer(); ++i) {
0174     layersPosZ_.push_back(rhtools_.getPositionLayer(i + 1).z());
0175     if (PatternRecognitionAlgoBaseT<TILES>::algo_verbosity_ > VerbosityLevel::Advanced) {
0176       edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "Layer " << i << " located at Z: " << layersPosZ_.back();
0177     }
0178   }
0179 
0180   clusters_.clear();
0181   tracksterSeedAlgoId_.clear();
0182 
0183   clusters_.resize(2 * rhtools_.lastLayer(false));
0184   std::vector<std::pair<int, int>> layerIdx2layerandSoa;  //used everywhere also to propagate cluster masking
0185 
0186   layerIdx2layerandSoa.reserve(input.layerClusters.size());
0187   unsigned int layerIdx = 0;
0188   for (auto const &lc : input.layerClusters) {
0189     if (input.mask[layerIdx] == 0.) {
0190       if (PatternRecognitionAlgoBaseT<TILES>::algo_verbosity_ > VerbosityLevel::Advanced) {
0191         edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "Skipping masked cluster: " << layerIdx;
0192       }
0193       layerIdx2layerandSoa.emplace_back(-1, -1);
0194       layerIdx++;
0195       continue;
0196     }
0197     const auto firstHitDetId = lc.hitsAndFractions()[0].first;
0198     int layer = rhtools_.getLayerWithOffset(firstHitDetId) - 1 +
0199                 rhtools_.lastLayer(false) * ((rhtools_.zside(firstHitDetId) + 1) >> 1);
0200     assert(layer >= 0);
0201     auto detId = lc.hitsAndFractions()[0].first;
0202     int layerClusterIndexInLayer = clusters_[layer].x.size();
0203     layerIdx2layerandSoa.emplace_back(layer, layerClusterIndexInLayer);
0204     float sum_x = 0.;
0205     float sum_y = 0.;
0206     float sum_sqr_x = 0.;
0207     float sum_sqr_y = 0.;
0208     float ref_x = lc.x();
0209     float ref_y = lc.y();
0210     float invClsize = 1. / lc.hitsAndFractions().size();
0211     for (auto const &hitsAndFractions : lc.hitsAndFractions()) {
0212       auto const &point = rhtools_.getPosition(hitsAndFractions.first);
0213       sum_x += point.x() - ref_x;
0214       sum_sqr_x += (point.x() - ref_x) * (point.x() - ref_x);
0215       sum_y += point.y() - ref_y;
0216       sum_sqr_y += (point.y() - ref_y) * (point.y() - ref_y);
0217     }
0218     // The variance of X for X uniform in circle of radius R is R^2/4,
0219     // therefore we multiply the sqrt(var) by 2 to have a rough estimate of the
0220     // radius. On the other hand, while averaging the x and y radius, we would
0221     // end up dividing by 2. Hence we omit the value here and in the average
0222     // below, too.
0223     float radius_x = sqrt((sum_sqr_x - (sum_x * sum_x) * invClsize) * invClsize);
0224     float radius_y = sqrt((sum_sqr_y - (sum_y * sum_y) * invClsize) * invClsize);
0225     if (PatternRecognitionAlgoBaseT<TILES>::algo_verbosity_ > VerbosityLevel::Advanced) {
0226       edm::LogVerbatim("PatternRecognitionbyCLUE3D")
0227           << "cluster rx: " << std::setw(5) << radius_x << ", ry: " << std::setw(5) << radius_y
0228           << ", r:  " << std::setw(5) << (radius_x + radius_y) << ", cells: " << std::setw(4)
0229           << lc.hitsAndFractions().size();
0230     }
0231 
0232     // The case of single cell layer clusters has to be handled differently.
0233 
0234     if (invClsize == 1.) {
0235       // Silicon case
0236       if (rhtools_.isSilicon(detId)) {
0237         radius_x = radius_y = rhtools_.getRadiusToSide(detId);
0238         if (PatternRecognitionAlgoBaseT<TILES>::algo_verbosity_ > VerbosityLevel::Advanced) {
0239           edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "Single cell cluster in silicon, rx: " << std::setw(5)
0240                                                          << radius_x << ", ry: " << std::setw(5) << radius_y;
0241         }
0242       } else {
0243         auto const &point = rhtools_.getPosition(detId);
0244         auto const &eta_phi_window = rhtools_.getScintDEtaDPhi(detId);
0245         radius_x = radius_y = point.perp() * eta_phi_window.second;
0246         if (PatternRecognitionAlgoBaseT<TILES>::algo_verbosity_ > VerbosityLevel::Advanced) {
0247           edm::LogVerbatim("PatternRecognitionbyCLUE3D")
0248               << "Single cell cluster in scintillator. rx: " << std::setw(5) << radius_x << ", ry: " << std::setw(5)
0249               << radius_y << ", eta-span: " << std::setw(5) << eta_phi_window.first << ", phi-span: " << std::setw(5)
0250               << eta_phi_window.second;
0251         }
0252       }
0253     }
0254 
0255     // Maybe check if these vectors can be reserved beforehands
0256     clusters_[layer].x.emplace_back(lc.x());
0257     clusters_[layer].y.emplace_back(lc.y());
0258     clusters_[layer].z.emplace_back(lc.z());
0259     clusters_[layer].r_over_absz.emplace_back(sqrt(lc.x() * lc.x() + lc.y() * lc.y()) / std::abs(lc.z()));
0260     clusters_[layer].radius.emplace_back(radius_x + radius_y);
0261     clusters_[layer].eta.emplace_back(lc.eta());
0262     clusters_[layer].phi.emplace_back(lc.phi());
0263     clusters_[layer].cells.push_back(lc.hitsAndFractions().size());
0264     clusters_[layer].algoId.push_back(lc.algo() - reco::CaloCluster::hgcal_em);
0265     clusters_[layer].isSilicon.push_back(rhtools_.isSilicon(detId));
0266     clusters_[layer].energy.emplace_back(lc.energy());
0267     clusters_[layer].isSeed.push_back(false);
0268     clusters_[layer].clusterIndex.emplace_back(-1);
0269     clusters_[layer].layerClusterOriginalIdx.emplace_back(layerIdx++);
0270     clusters_[layer].nearestHigher.emplace_back(-1, -1);
0271     clusters_[layer].rho.emplace_back(0.f);
0272     clusters_[layer].z_extension.emplace_back(0.f);
0273     clusters_[layer].delta.emplace_back(
0274         std::make_pair(std::numeric_limits<float>::max(), std::numeric_limits<int>::max()));
0275   }
0276   for (unsigned int layer = 0; layer < clusters_.size(); layer++) {
0277     clusters_[layer].followers.resize(clusters_[layer].x.size());
0278   }
0279 
0280   auto lastLayerPerSide = static_cast<int>(rhtools_.lastLayer(false));
0281   int maxLayer = 2 * lastLayerPerSide - 1;
0282   std::vector<int> numberOfClustersPerLayer(maxLayer, 0);
0283   for (int i = 0; i <= maxLayer; i++) {
0284     calculateLocalDensity(input.tiles, i, layerIdx2layerandSoa);
0285   }
0286   for (int i = 0; i <= maxLayer; i++) {
0287     calculateDistanceToHigher(input.tiles, i, layerIdx2layerandSoa);
0288   }
0289 
0290   auto nTracksters = findAndAssignTracksters(input.tiles, layerIdx2layerandSoa);
0291   if (PatternRecognitionAlgoBaseT<TILES>::algo_verbosity_ > VerbosityLevel::Advanced) {
0292     edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "Reconstructed " << nTracksters << " tracksters" << std::endl;
0293     dumpClusters(input.tiles, layerIdx2layerandSoa, eventNumber);
0294   }
0295 
0296   // Build Trackster
0297   result.resize(nTracksters);
0298 
0299   for (unsigned int layer = 0; layer < clusters_.size(); ++layer) {
0300     const auto &thisLayer = clusters_[layer];
0301     if (PatternRecognitionAlgoBaseT<TILES>::algo_verbosity_ > VerbosityLevel::Advanced) {
0302       edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "Examining Layer: " << layer;
0303     }
0304     for (unsigned int lc = 0; lc < thisLayer.x.size(); ++lc) {
0305       if (PatternRecognitionAlgoBaseT<TILES>::algo_verbosity_ > VerbosityLevel::Advanced) {
0306         edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "Trackster " << thisLayer.clusterIndex[lc];
0307       }
0308       if (thisLayer.clusterIndex[lc] >= 0) {
0309         if (PatternRecognitionAlgoBaseT<TILES>::algo_verbosity_ > VerbosityLevel::Advanced) {
0310           edm::LogVerbatim("PatternRecognitionbyCLUE3D") << " adding lcIdx: " << thisLayer.layerClusterOriginalIdx[lc];
0311         }
0312         result[thisLayer.clusterIndex[lc]].vertices().push_back(thisLayer.layerClusterOriginalIdx[lc]);
0313         result[thisLayer.clusterIndex[lc]].vertex_multiplicity().push_back(1);
0314         // loop over followers
0315         for (auto [follower_lyrIdx, follower_soaIdx] : thisLayer.followers[lc]) {
0316           std::array<unsigned int, 2> edge = {
0317               {(unsigned int)thisLayer.layerClusterOriginalIdx[lc],
0318                (unsigned int)clusters_[follower_lyrIdx].layerClusterOriginalIdx[follower_soaIdx]}};
0319           result[thisLayer.clusterIndex[lc]].edges().push_back(edge);
0320         }
0321       }
0322     }
0323   }
0324   size_t tracksterIndex = 0;
0325   result.erase(std::remove_if(std::begin(result),
0326                               std::end(result),
0327                               [&](auto const &v) {
0328                                 return static_cast<int>(v.vertices().size()) <
0329                                        minNumLayerCluster_[tracksterSeedAlgoId_[tracksterIndex++]];
0330                               }),
0331                result.end());
0332   if (doPidCut_) {
0333     energyRegressionAndID(input.layerClusters, input.tfSession, result);
0334     result.erase(std::remove_if(std::begin(result),
0335                                 std::end(result),
0336                                 [&](auto const &v) {
0337                                   auto const &hadProb =
0338                                       v.id_probability(ticl::Trackster::ParticleType::charged_hadron) +
0339                                       v.id_probability(ticl::Trackster::ParticleType::neutral_hadron);
0340                                   return hadProb >= cutHadProb_;
0341                                 }),
0342                  result.end());
0343   }
0344   result.shrink_to_fit();
0345 
0346   ticl::assignPCAtoTracksters(result,
0347                               input.layerClusters,
0348                               input.layerClustersTime,
0349                               rhtools_.getPositionLayer(rhtools_.lastLayerEE(false), false).z(),
0350                               computeLocalTime_);
0351 
0352   if (PatternRecognitionAlgoBaseT<TILES>::algo_verbosity_ > VerbosityLevel::Advanced) {
0353     for (auto const &t : result) {
0354       edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "Barycenter: " << t.barycenter();
0355       edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "LCs: " << t.vertices().size();
0356       edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "Energy: " << t.raw_energy();
0357       edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "Regressed: " << t.regressed_energy();
0358     }
0359   }
0360 
0361   // Dump Tracksters information
0362   if (PatternRecognitionAlgoBaseT<TILES>::algo_verbosity_ > VerbosityLevel::Advanced) {
0363     dumpTracksters(layerIdx2layerandSoa, eventNumber, result);
0364   }
0365 
0366   // Reset internal clusters_ structure of array for next event
0367   reset();
0368   if (PatternRecognitionAlgoBaseT<TILES>::algo_verbosity_ > VerbosityLevel::Advanced) {
0369     edm::LogVerbatim("PatternRecognitionbyCLUE3D") << std::endl;
0370   }
0371 }
0372 
0373 template <typename TILES>
0374 void PatternRecognitionbyCLUE3D<TILES>::energyRegressionAndID(const std::vector<reco::CaloCluster> &layerClusters,
0375                                                               const tensorflow::Session *eidSession,
0376                                                               std::vector<Trackster> &tracksters) {
0377   // Energy regression and particle identification strategy:
0378   //
0379   // 1. Set default values for regressed energy and particle id for each trackster.
0380   // 2. Store indices of tracksters whose total sum of cluster energies is above the
0381   //    eidMinClusterEnergy_ (GeV) threshold. Inference is not applied for soft tracksters.
0382   // 3. When no trackster passes the selection, return.
0383   // 4. Create input and output tensors. The batch dimension is determined by the number of
0384   //    selected tracksters.
0385   // 5. Fill input tensors with layer cluster features. Per layer, clusters are ordered descending
0386   //    by energy. Given that tensor data is contiguous in memory, we can use pointer arithmetic to
0387   //    fill values, even with batching.
0388   // 6. Zero-fill features for empty clusters in each layer.
0389   // 7. Batched inference.
0390   // 8. Assign the regressed energy and id probabilities to each trackster.
0391   //
0392   // Indices used throughout this method:
0393   // i -> batch element / trackster
0394   // j -> layer
0395   // k -> cluster
0396   // l -> feature
0397 
0398   // set default values per trackster, determine if the cluster energy threshold is passed,
0399   // and store indices of hard tracksters
0400   std::vector<int> tracksterIndices;
0401   for (int i = 0; i < static_cast<int>(tracksters.size()); i++) {
0402     // calculate the cluster energy sum (2)
0403     // note: after the loop, sumClusterEnergy might be just above the threshold which is enough to
0404     // decide whether to run inference for the trackster or not
0405     float sumClusterEnergy = 0.;
0406     for (const unsigned int &vertex : tracksters[i].vertices()) {
0407       sumClusterEnergy += static_cast<float>(layerClusters[vertex].energy());
0408       // there might be many clusters, so try to stop early
0409       if (sumClusterEnergy >= eidMinClusterEnergy_) {
0410         // set default values (1)
0411         tracksters[i].setRegressedEnergy(0.f);
0412         tracksters[i].zeroProbabilities();
0413         tracksterIndices.push_back(i);
0414         break;
0415       }
0416     }
0417   }
0418 
0419   // do nothing when no trackster passes the selection (3)
0420   int batchSize = static_cast<int>(tracksterIndices.size());
0421   if (batchSize == 0) {
0422     return;
0423   }
0424 
0425   // create input and output tensors (4)
0426   tensorflow::TensorShape shape({batchSize, eidNLayers_, eidNClusters_, eidNFeatures_});
0427   tensorflow::Tensor input(tensorflow::DT_FLOAT, shape);
0428   tensorflow::NamedTensorList inputList = {{eidInputName_, input}};
0429 
0430   std::vector<tensorflow::Tensor> outputs;
0431   std::vector<std::string> outputNames;
0432   if (!eidOutputNameEnergy_.empty()) {
0433     outputNames.push_back(eidOutputNameEnergy_);
0434   }
0435   if (!eidOutputNameId_.empty()) {
0436     outputNames.push_back(eidOutputNameId_);
0437   }
0438 
0439   // fill input tensor (5)
0440   for (int i = 0; i < batchSize; i++) {
0441     const Trackster &trackster = tracksters[tracksterIndices[i]];
0442 
0443     // per layer, we only consider the first eidNClusters_ clusters in terms of energy, so in order
0444     // to avoid creating large / nested structures to do the sorting for an unknown number of total
0445     // clusters, create a sorted list of layer cluster indices to keep track of the filled clusters
0446     std::vector<int> clusterIndices(trackster.vertices().size());
0447     for (int k = 0; k < (int)trackster.vertices().size(); k++) {
0448       clusterIndices[k] = k;
0449     }
0450     sort(clusterIndices.begin(), clusterIndices.end(), [&layerClusters, &trackster](const int &a, const int &b) {
0451       return layerClusters[trackster.vertices(a)].energy() > layerClusters[trackster.vertices(b)].energy();
0452     });
0453 
0454     // keep track of the number of seen clusters per layer
0455     std::vector<int> seenClusters(eidNLayers_);
0456 
0457     // loop through clusters by descending energy
0458     for (const int &k : clusterIndices) {
0459       // get features per layer and cluster and store the values directly in the input tensor
0460       const reco::CaloCluster &cluster = layerClusters[trackster.vertices(k)];
0461       int j = rhtools_.getLayerWithOffset(cluster.hitsAndFractions()[0].first) - 1;
0462       if (j < eidNLayers_ && seenClusters[j] < eidNClusters_) {
0463         // get the pointer to the first feature value for the current batch, layer and cluster
0464         float *features = &input.tensor<float, 4>()(i, j, seenClusters[j], 0);
0465 
0466         // fill features
0467         *(features++) = float(cluster.energy() / float(trackster.vertex_multiplicity(k)));
0468         *(features++) = float(std::abs(cluster.eta()));
0469         *(features) = float(cluster.phi());
0470 
0471         // increment seen clusters
0472         seenClusters[j]++;
0473       }
0474     }
0475 
0476     // zero-fill features of empty clusters in each layer (6)
0477     for (int j = 0; j < eidNLayers_; j++) {
0478       for (int k = seenClusters[j]; k < eidNClusters_; k++) {
0479         float *features = &input.tensor<float, 4>()(i, j, k, 0);
0480         for (int l = 0; l < eidNFeatures_; l++) {
0481           *(features++) = 0.f;
0482         }
0483       }
0484     }
0485   }
0486 
0487   // run the inference (7)
0488   tensorflow::run(eidSession, inputList, outputNames, &outputs);
0489 
0490   // store regressed energy per trackster (8)
0491   if (!eidOutputNameEnergy_.empty()) {
0492     // get the pointer to the energy tensor, dimension is batch x 1
0493     float *energy = outputs[0].flat<float>().data();
0494 
0495     for (const int &i : tracksterIndices) {
0496       tracksters[i].setRegressedEnergy(*(energy++));
0497     }
0498   }
0499 
0500   // store id probabilities per trackster (8)
0501   if (!eidOutputNameId_.empty()) {
0502     // get the pointer to the id probability tensor, dimension is batch x id_probabilities.size()
0503     int probsIdx = eidOutputNameEnergy_.empty() ? 0 : 1;
0504     float *probs = outputs[probsIdx].flat<float>().data();
0505 
0506     for (const int &i : tracksterIndices) {
0507       tracksters[i].setProbabilities(probs);
0508       probs += tracksters[i].id_probabilities().size();
0509     }
0510   }
0511 }
0512 
0513 template <typename TILES>
0514 void PatternRecognitionbyCLUE3D<TILES>::calculateLocalDensity(
0515     const TILES &tiles, const int layerId, const std::vector<std::pair<int, int>> &layerIdx2layerandSoa) {
0516   constexpr int nEtaBin = TILES::constants_type_t::nEtaBins;
0517   constexpr int nPhiBin = TILES::constants_type_t::nPhiBins;
0518   auto &clustersOnLayer = clusters_[layerId];
0519   unsigned int numberOfClusters = clustersOnLayer.x.size();
0520 
0521   auto isReachable = [](float r0, float r1, float phi0, float phi1, float delta_sqr) -> bool {
0522     auto delta_phi = reco::deltaPhi(phi0, phi1);
0523     return (r0 - r1) * (r0 - r1) + r1 * r1 * delta_phi * delta_phi < delta_sqr;
0524   };
0525   auto distance_debug = [&](float x1, float x2, float y1, float y2) -> float {
0526     return sqrt((x1 - x2) * (x1 - x2) + (y1 - y2) * (y1 - y2));
0527   };
0528 
0529   for (unsigned int i = 0; i < numberOfClusters; i++) {
0530     auto algoId = clustersOnLayer.algoId[i];
0531     // We need to partition the two sides of the HGCAL detector
0532     auto lastLayerPerSide = static_cast<int>(rhtools_.lastLayer(false));
0533     int minLayer = 0;
0534     int maxLayer = 2 * lastLayerPerSide - 1;
0535     if (layerId < lastLayerPerSide) {
0536       minLayer = std::max(layerId - densitySiblingLayers_[algoId], minLayer);
0537       maxLayer = std::min(layerId + densitySiblingLayers_[algoId], lastLayerPerSide - 1);
0538     } else {
0539       minLayer = std::max(layerId - densitySiblingLayers_[algoId], lastLayerPerSide);
0540       maxLayer = std::min(layerId + densitySiblingLayers_[algoId], maxLayer);
0541     }
0542     float deltaLayersZ = std::abs(layersPosZ_[maxLayer % lastLayerPerSide] - layersPosZ_[minLayer % lastLayerPerSide]);
0543 
0544     for (int currentLayer = minLayer; currentLayer <= maxLayer; currentLayer++) {
0545       if (PatternRecognitionAlgoBaseT<TILES>::algo_verbosity_ > VerbosityLevel::Advanced) {
0546         edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "RefLayer: " << layerId << " SoaIDX: " << i;
0547         edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "NextLayer: " << currentLayer;
0548       }
0549       const auto &tileOnLayer = tiles[currentLayer];
0550       bool onSameLayer = (currentLayer == layerId);
0551       if (PatternRecognitionAlgoBaseT<TILES>::algo_verbosity_ > VerbosityLevel::Advanced) {
0552         edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "onSameLayer: " << onSameLayer;
0553       }
0554       const int etaWindow = 2;
0555       const int phiWindow = 2;
0556       int etaBinMin = std::max(tileOnLayer.etaBin(clustersOnLayer.eta[i]) - etaWindow, 0);
0557       int etaBinMax = std::min(tileOnLayer.etaBin(clustersOnLayer.eta[i]) + etaWindow, nEtaBin - 1);
0558       int phiBinMin = tileOnLayer.phiBin(clustersOnLayer.phi[i]) - phiWindow;
0559       int phiBinMax = tileOnLayer.phiBin(clustersOnLayer.phi[i]) + phiWindow;
0560       if (PatternRecognitionAlgoBaseT<TILES>::algo_verbosity_ > VerbosityLevel::Advanced) {
0561         edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "eta: " << clustersOnLayer.eta[i];
0562         edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "phi: " << clustersOnLayer.phi[i];
0563         edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "etaBinMin: " << etaBinMin << ", etaBinMax: " << etaBinMax;
0564         edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "phiBinMin: " << phiBinMin << ", phiBinMax: " << phiBinMax;
0565       }
0566       for (int ieta = etaBinMin; ieta <= etaBinMax; ++ieta) {
0567         auto offset = ieta * nPhiBin;
0568         if (PatternRecognitionAlgoBaseT<TILES>::algo_verbosity_ > VerbosityLevel::Advanced) {
0569           edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "offset: " << offset;
0570         }
0571         for (int iphi_it = phiBinMin; iphi_it <= phiBinMax; ++iphi_it) {
0572           int iphi = ((iphi_it % nPhiBin + nPhiBin) % nPhiBin);
0573           if (PatternRecognitionAlgoBaseT<TILES>::algo_verbosity_ > VerbosityLevel::Advanced) {
0574             edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "iphi: " << iphi;
0575             edm::LogVerbatim("PatternRecognitionbyCLUE3D")
0576                 << "Entries in tileBin: " << tileOnLayer[offset + iphi].size();
0577           }
0578           for (auto otherClusterIdx : tileOnLayer[offset + iphi]) {
0579             auto const &layerandSoa = layerIdx2layerandSoa[otherClusterIdx];
0580             // Skip masked layer clusters
0581             if ((layerandSoa.first == -1) && (layerandSoa.second == -1)) {
0582               if (PatternRecognitionAlgoBaseT<TILES>::algo_verbosity_ > VerbosityLevel::Advanced) {
0583                 edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "Skipping masked layerIdx " << otherClusterIdx;
0584               }
0585               continue;
0586             }
0587             auto const &clustersLayer = clusters_[layerandSoa.first];
0588             if (PatternRecognitionAlgoBaseT<TILES>::algo_verbosity_ > VerbosityLevel::Advanced) {
0589               edm::LogVerbatim("PatternRecognitionbyCLUE3D")
0590                   << "OtherLayer: " << layerandSoa.first << " SoaIDX: " << layerandSoa.second;
0591               edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "OtherEta: " << clustersLayer.eta[layerandSoa.second];
0592               edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "OtherPhi: " << clustersLayer.phi[layerandSoa.second];
0593             }
0594 
0595             bool onSameCluster = clustersOnLayer.layerClusterOriginalIdx[i] == otherClusterIdx;
0596             if (onSameLayer && !densityOnSameLayer_ && !onSameCluster) {
0597               if (PatternRecognitionAlgoBaseT<TILES>::algo_verbosity_ > VerbosityLevel::Advanced) {
0598                 edm::LogVerbatim("PatternRecognitionbyCLUE3D")
0599                     << "Skipping different cluster " << otherClusterIdx << "in the same layer " << currentLayer;
0600               }
0601               continue;
0602             }
0603 
0604             bool reachable = false;
0605             if (useAbsoluteProjectiveScale_) {
0606               if (useClusterDimensionXY_) {
0607                 reachable = isReachable(clustersOnLayer.r_over_absz[i] * clustersOnLayer.z[i],
0608                                         clustersLayer.r_over_absz[layerandSoa.second] * clustersOnLayer.z[i],
0609                                         clustersOnLayer.phi[i],
0610                                         clustersLayer.phi[layerandSoa.second],
0611                                         clustersOnLayer.radius[i] * clustersOnLayer.radius[i]);
0612               } else {
0613                 // Still differentiate between silicon and Scintillator.
0614                 // Scintillator has yet to be studied further.
0615                 if (clustersOnLayer.isSilicon[i]) {
0616                   reachable = isReachable(clustersOnLayer.r_over_absz[i] * clustersOnLayer.z[i],
0617                                           clustersLayer.r_over_absz[layerandSoa.second] * clustersOnLayer.z[i],
0618                                           clustersOnLayer.phi[i],
0619                                           clustersLayer.phi[layerandSoa.second],
0620                                           densityXYDistanceSqr_[algoId]);
0621                 } else {
0622                   reachable = isReachable(clustersOnLayer.r_over_absz[i] * clustersOnLayer.z[i],
0623                                           clustersLayer.r_over_absz[layerandSoa.second] * clustersOnLayer.z[i],
0624                                           clustersOnLayer.phi[i],
0625                                           clustersLayer.phi[layerandSoa.second],
0626                                           clustersOnLayer.radius[i] * clustersOnLayer.radius[i]);
0627                 }
0628               }
0629             } else {
0630               reachable = (reco::deltaR2(clustersOnLayer.eta[i],
0631                                          clustersOnLayer.phi[i],
0632                                          clustersLayer.eta[layerandSoa.second],
0633                                          clustersLayer.phi[layerandSoa.second]) < densityEtaPhiDistanceSqr_[algoId]);
0634             }
0635             if (PatternRecognitionAlgoBaseT<TILES>::algo_verbosity_ > VerbosityLevel::Advanced) {
0636               edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "Distance[eta,phi]: "
0637                                                              << reco::deltaR2(clustersOnLayer.eta[i],
0638                                                                               clustersOnLayer.phi[i],
0639                                                                               clustersLayer.eta[layerandSoa.second],
0640                                                                               clustersLayer.phi[layerandSoa.second]);
0641               auto dist = distance_debug(
0642                   clustersOnLayer.r_over_absz[i],
0643                   clustersLayer.r_over_absz[layerandSoa.second],
0644                   clustersOnLayer.r_over_absz[i] * std::abs(clustersOnLayer.phi[i]),
0645                   clustersLayer.r_over_absz[layerandSoa.second] * std::abs(clustersLayer.phi[layerandSoa.second]));
0646               edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "Distance[cm]: " << (dist * clustersOnLayer.z[i]);
0647               edm::LogVerbatim("PatternRecognitionbyCLUE3D")
0648                   << "Energy Other:   " << clustersLayer.energy[layerandSoa.second];
0649               edm::LogVerbatim("PatternRecognitionbyCLUE3D") << "Cluster radius: " << clustersOnLayer.radius[i];
0650             }
0651             if (reachable) {
0652               float factor_same_layer_different_cluster = (onSameLayer && !densityOnSameLayer_) ? 0.f : 1.f;
0653               auto energyToAdd = (clustersOnLayer.layerClusterOriginalIdx[i] == otherClusterIdx
0654                                       ? 1.f
0655                                       : kernelDensityFactor_[algoId] * factor_same_layer_different_cluster) *
0656                                  clustersLayer.energy[layerandSoa.second];
0657               clustersOnLayer.rho[i] += energyToAdd;
0658               clustersOnLayer.z_extension[i] = deltaLayersZ;
0659               if (PatternRecognitionAlgoBaseT<TILES>::algo_verbosity_ > VerbosityLevel::Advanced) {
0660                 edm::LogVerbatim("PatternRecognitionbyCLUE3D")
0661                     << "Adding " << energyToAdd << " partial " << clustersOnLayer.rho[i];
0662               }
0663             }
0664           }  // end of loop on possible compatible clusters
0665         }    // end of loop over phi-bin region
0666       }      // end of loop over eta-bin region
0667     }        // end of loop on the sibling layers
0668     if (rescaleDensityByZ_) {
0669       if (PatternRecognitionAlgoBaseT<TILES>::algo_verbosity_ > VerbosityLevel::Advanced) {
0670         edm::LogVerbatim("PatternRecognitionbyCLUE3D")
0671             << "Rescaling original density: " << clustersOnLayer.rho[i] << " by Z: " << deltaLayersZ
0672             << " to final density/cm: " << clustersOnLayer.rho[i] / deltaLayersZ;
0673       }
0674       clustersOnLayer.rho[i] /= deltaLayersZ;
0675     }
0676   }  // end of loop over clusters on this layer
0677 }
0678 
0679 template <typename TILES>
0680 void PatternRecognitionbyCLUE3D<TILES>::calculateDistanceToHigher(
0681     const TILES &tiles, const int layerId, const std::vector<std::pair<int, int>> &layerIdx2layerandSoa) {
0682   constexpr int nEtaBin = TILES::constants_type_t::nEtaBins;
0683   constexpr int nPhiBin = TILES::constants_type_t::nPhiBins;
0684   auto &clustersOnLayer = clusters_[layerId];
0685   unsigned int numberOfClusters = clustersOnLayer.x.size();
0686 
0687   auto distanceSqr = [](float r0, float r1, float phi0, float phi1) -> float {
0688     auto delta_phi = reco::deltaPhi(phi0, phi1);
0689     return (r0 - r1) * (r0 - r1) + r1 * r1 * delta_phi * delta_phi;
0690   };
0691 
0692   for (unsigned int i = 0; i < numberOfClusters; i++) {
0693     if (PatternRecognitionAlgoBaseT<TILES>::algo_verbosity_ > VerbosityLevel::Advanced) {
0694       edm::LogVerbatim("PatternRecognitionbyCLUE3D")
0695           << "Starting searching nearestHigher on " << layerId << " with rho: " << clustersOnLayer.rho[i]
0696           << " at eta, phi: " << tiles[layerId].etaBin(clustersOnLayer.eta[i]) << ", "
0697           << tiles[layerId].phiBin(clustersOnLayer.phi[i]);
0698     }
0699     // We need to partition the two sides of the HGCAL detector
0700     auto lastLayerPerSide = static_cast<int>(rhtools_.lastLayer(false));
0701     int minLayer = 0;
0702     int maxLayer = 2 * lastLayerPerSide - 1;
0703     auto algoId = clustersOnLayer.algoId[i];
0704     if (layerId < lastLayerPerSide) {
0705       minLayer = std::max(layerId - densitySiblingLayers_[algoId], minLayer);
0706       maxLayer = std::min(layerId + densitySiblingLayers_[algoId], lastLayerPerSide - 1);
0707     } else {
0708       minLayer = std::max(layerId - densitySiblingLayers_[algoId], lastLayerPerSide + 1);
0709       maxLayer = std::min(layerId + densitySiblingLayers_[algoId], maxLayer);
0710     }
0711     constexpr float maxDelta = std::numeric_limits<float>::max();
0712     float i_delta = maxDelta;
0713     std::pair<int, int> i_nearestHigher(-1, -1);
0714     std::pair<float, int> nearest_distances(maxDelta, std::numeric_limits<int>::max());
0715     for (int currentLayer = minLayer; currentLayer <= maxLayer; currentLayer++) {
0716       if (!nearestHigherOnSameLayer_ && (layerId == currentLayer))
0717         continue;
0718       const auto &tileOnLayer = tiles[currentLayer];
0719       int etaWindow = 1;
0720       int phiWindow = 1;
0721       int etaBinMin = std::max(tileOnLayer.etaBin(clustersOnLayer.eta[i]) - etaWindow, 0);
0722       int etaBinMax = std::min(tileOnLayer.etaBin(clustersOnLayer.eta[i]) + etaWindow, nEtaBin - 1);
0723       int phiBinMin = tileOnLayer.phiBin(clustersOnLayer.phi[i]) - phiWindow;
0724       int phiBinMax = tileOnLayer.phiBin(clustersOnLayer.phi[i]) + phiWindow;
0725       for (int ieta = etaBinMin; ieta <= etaBinMax; ++ieta) {
0726         auto offset = ieta * nPhiBin;
0727         for (int iphi_it = phiBinMin; iphi_it <= phiBinMax; ++iphi_it) {
0728           int iphi = ((iphi_it % nPhiBin + nPhiBin) % nPhiBin);
0729           if (PatternRecognitionAlgoBaseT<TILES>::algo_verbosity_ > VerbosityLevel::Advanced) {
0730             edm::LogVerbatim("PatternRecognitionbyCLUE3D")
0731                 << "Searching nearestHigher on " << currentLayer << " eta, phi: " << ieta << ", " << iphi_it << " "
0732                 << iphi << " " << offset << " " << (offset + iphi);
0733           }
0734           for (auto otherClusterIdx : tileOnLayer[offset + iphi]) {
0735             auto const &layerandSoa = layerIdx2layerandSoa[otherClusterIdx];
0736             // Skip masked layer clusters
0737             if ((layerandSoa.first == -1) && (layerandSoa.second == -1))
0738               continue;
0739             auto const &clustersOnOtherLayer = clusters_[layerandSoa.first];
0740             auto dist = maxDelta;
0741             auto dist_transverse = maxDelta;
0742             int dist_layers = std::abs(layerandSoa.first - layerId);
0743             if (useAbsoluteProjectiveScale_) {
0744               dist_transverse = distanceSqr(clustersOnLayer.r_over_absz[i] * clustersOnLayer.z[i],
0745                                             clustersOnOtherLayer.r_over_absz[layerandSoa.second] * clustersOnLayer.z[i],
0746                                             clustersOnLayer.phi[i],
0747                                             clustersOnOtherLayer.phi[layerandSoa.second]);
0748               // Add Z-scale to the final distance
0749               dist = dist_transverse;
0750             } else {
0751               dist = reco::deltaR2(clustersOnLayer.eta[i],
0752                                    clustersOnLayer.phi[i],
0753                                    clustersOnOtherLayer.eta[layerandSoa.second],
0754                                    clustersOnOtherLayer.phi[layerandSoa.second]);
0755               dist_transverse = dist;
0756             }
0757             bool foundHigher = (clustersOnOtherLayer.rho[layerandSoa.second] > clustersOnLayer.rho[i]) ||
0758                                (clustersOnOtherLayer.rho[layerandSoa.second] == clustersOnLayer.rho[i] &&
0759                                 clustersOnOtherLayer.layerClusterOriginalIdx[layerandSoa.second] >
0760                                     clustersOnLayer.layerClusterOriginalIdx[i]);
0761             if (PatternRecognitionAlgoBaseT<TILES>::algo_verbosity_ > VerbosityLevel::Advanced) {
0762               edm::LogVerbatim("PatternRecognitionbyCLUE3D")
0763                   << "Searching nearestHigher on " << currentLayer
0764                   << " with rho: " << clustersOnOtherLayer.rho[layerandSoa.second]
0765                   << " on layerIdxInSOA: " << layerandSoa.first << ", " << layerandSoa.second
0766                   << " with distance: " << sqrt(dist) << " foundHigher: " << foundHigher;
0767             }
0768             if (foundHigher && dist <= i_delta) {
0769               // update i_delta
0770               i_delta = dist;
0771               nearest_distances = std::make_pair(sqrt(dist_transverse), dist_layers);
0772               // update i_nearestHigher
0773               i_nearestHigher = layerandSoa;
0774             }
0775           }  // End of loop on clusters
0776         }    // End of loop on phi bins
0777       }      // End of loop on eta bins
0778     }        // End of loop on layers
0779 
0780     clustersOnLayer.delta[i] = nearest_distances;
0781     clustersOnLayer.nearestHigher[i] = i_nearestHigher;
0782   }  // End of loop on clusters
0783 }
0784 
0785 template <typename TILES>
0786 int PatternRecognitionbyCLUE3D<TILES>::findAndAssignTracksters(
0787     const TILES &tiles, const std::vector<std::pair<int, int>> &layerIdx2layerandSoa) {
0788   unsigned int nTracksters = 0;
0789   std::vector<std::pair<int, int>> localStack;
0790   const auto &critical_transverse_distance =
0791       useAbsoluteProjectiveScale_ ? criticalXYDistance_ : criticalEtaPhiDistance_;
0792   // find cluster seeds and outlier
0793   for (unsigned int layer = 0; layer < 2 * rhtools_.lastLayer(); layer++) {
0794     auto &clustersOnLayer = clusters_[layer];
0795     unsigned int numberOfClusters = clustersOnLayer.x.size();
0796     for (unsigned int i = 0; i < numberOfClusters; i++) {
0797       // initialize clusterIndex
0798       clustersOnLayer.clusterIndex[i] = -1;
0799       auto algoId = clustersOnLayer.algoId[i];
0800       bool isSeed = (clustersOnLayer.delta[i].first > critical_transverse_distance[algoId] ||
0801                      clustersOnLayer.delta[i].second > criticalZDistanceLyr_[algoId]) &&
0802                     (clustersOnLayer.rho[i] >= criticalDensity_[algoId]) &&
0803                     (clustersOnLayer.energy[i] / clustersOnLayer.rho[i] > criticalSelfDensity_[algoId]);
0804       if (!clustersOnLayer.isSilicon[i]) {
0805         isSeed = (clustersOnLayer.delta[i].first > clustersOnLayer.radius[i] ||
0806                   clustersOnLayer.delta[i].second > criticalZDistanceLyr_[algoId]) &&
0807                  (clustersOnLayer.rho[i] >= criticalDensity_[algoId]) &&
0808                  (clustersOnLayer.energy[i] / clustersOnLayer.rho[i] > criticalSelfDensity_[algoId]);
0809       }
0810       bool isOutlier =
0811           (clustersOnLayer.delta[i].first > outlierMultiplier_[algoId] * critical_transverse_distance[algoId]) &&
0812           (clustersOnLayer.rho[i] < criticalDensity_[algoId]);
0813       if (isSeed) {
0814         if (PatternRecognitionAlgoBaseT<TILES>::algo_verbosity_ > VerbosityLevel::Advanced) {
0815           edm::LogVerbatim("PatternRecognitionbyCLUE3D")
0816               << "Found seed on Layer " << layer << " SOAidx: " << i << " assigned ClusterIdx: " << nTracksters;
0817         }
0818         clustersOnLayer.clusterIndex[i] = nTracksters++;
0819         tracksterSeedAlgoId_.push_back(algoId);
0820         clustersOnLayer.isSeed[i] = true;
0821         localStack.emplace_back(layer, i);
0822       } else if (!isOutlier) {
0823         auto [lyrIdx, soaIdx] = clustersOnLayer.nearestHigher[i];
0824         if (PatternRecognitionAlgoBaseT<TILES>::algo_verbosity_ > VerbosityLevel::Advanced) {
0825           edm::LogVerbatim("PatternRecognitionbyCLUE3D")
0826               << "Found follower on Layer " << layer << " SOAidx: " << i << " attached to cluster on layer: " << lyrIdx
0827               << " SOAidx: " << soaIdx;
0828         }
0829         if (lyrIdx >= 0)
0830           clusters_[lyrIdx].followers[soaIdx].emplace_back(layer, i);
0831       } else {
0832         if (PatternRecognitionAlgoBaseT<TILES>::algo_verbosity_ > VerbosityLevel::Advanced) {
0833           edm::LogVerbatim("PatternRecognitionbyCLUE3D")
0834               << "Found Outlier on Layer " << layer << " SOAidx: " << i << " with rho: " << clustersOnLayer.rho[i]
0835               << " and delta: " << clustersOnLayer.delta[i].first << ", " << clustersOnLayer.delta[i].second;
0836         }
0837       }
0838     }
0839   }
0840 
0841   // Propagate cluster index
0842   while (!localStack.empty()) {
0843     auto [lyrIdx, soaIdx] = localStack.back();
0844     auto &thisSeed = clusters_[lyrIdx].followers[soaIdx];
0845     localStack.pop_back();
0846 
0847     // loop over followers
0848     for (auto [follower_lyrIdx, follower_soaIdx] : thisSeed) {
0849       // pass id to a follower
0850       clusters_[follower_lyrIdx].clusterIndex[follower_soaIdx] = clusters_[lyrIdx].clusterIndex[soaIdx];
0851       // push this follower to localStack
0852       localStack.emplace_back(follower_lyrIdx, follower_soaIdx);
0853     }
0854   }
0855   return nTracksters;
0856 }
0857 
0858 template <typename TILES>
0859 void PatternRecognitionbyCLUE3D<TILES>::fillPSetDescription(edm::ParameterSetDescription &iDesc) {
0860   iDesc.add<int>("algo_verbosity", 0);
0861   iDesc.add<std::vector<double>>("criticalDensity", {4, 4, 4})->setComment("in GeV");
0862   iDesc.add<std::vector<double>>("criticalSelfDensity", {0.15, 0.15, 0.15} /* roughly 1/(densitySiblingLayers+1) */)
0863       ->setComment("Minimum ratio of self_energy/local_density to become a seed.");
0864   iDesc.add<std::vector<int>>("densitySiblingLayers", {3, 3, 3})
0865       ->setComment(
0866           "inclusive, layers to consider while computing local density and searching for nearestHigher higher");
0867   iDesc.add<std::vector<double>>("densityEtaPhiDistanceSqr", {0.0008, 0.0008, 0.0008})
0868       ->setComment("in eta,phi space, distance to consider for local density");
0869   iDesc.add<std::vector<double>>("densityXYDistanceSqr", {3.24, 3.24, 3.24})
0870       ->setComment("in cm, distance on the transverse plane to consider for local density");
0871   iDesc.add<std::vector<double>>("kernelDensityFactor", {0.2, 0.2, 0.2})
0872       ->setComment("Kernel factor to be applied to other LC while computing the local density");
0873   iDesc.add<bool>("densityOnSameLayer", false);
0874   iDesc.add<bool>("nearestHigherOnSameLayer", false)
0875       ->setComment("Allow the nearestHigher to be located on the same layer");
0876   iDesc.add<bool>("useAbsoluteProjectiveScale", true)
0877       ->setComment("Express all cuts in terms of r/z*z_0{,phi} projective variables");
0878   iDesc.add<bool>("useClusterDimensionXY", false)
0879       ->setComment(
0880           "Boolean. If true use the estimated cluster radius to determine the cluster compatibility while computing "
0881           "the local density");
0882   iDesc.add<bool>("rescaleDensityByZ", false)
0883       ->setComment(
0884           "Rescale local density by the extension of the Z 'volume' explored. The transvere dimension is, at present, "
0885           "fixed and factored out.");
0886   iDesc.add<std::vector<double>>("criticalEtaPhiDistance", {0.025, 0.025, 0.025})
0887       ->setComment("Minimal distance in eta,phi space from nearestHigher to become a seed");
0888   iDesc.add<std::vector<double>>("criticalXYDistance", {1.8, 1.8, 1.8})
0889       ->setComment("Minimal distance in cm on the XY plane from nearestHigher to become a seed");
0890   iDesc.add<std::vector<int>>("criticalZDistanceLyr", {5, 5, 5})
0891       ->setComment("Minimal distance in layers along the Z axis from nearestHigher to become a seed");
0892   iDesc.add<std::vector<double>>("outlierMultiplier", {2, 2, 2})
0893       ->setComment("Minimal distance in transverse space from nearestHigher to become an outlier");
0894   iDesc.add<std::vector<int>>("minNumLayerCluster", {2, 2, 2})->setComment("Not Inclusive");
0895   iDesc.add<bool>("doPidCut", false);
0896   iDesc.add<double>("cutHadProb", 0.5);
0897   iDesc.add<std::string>("eid_input_name", "input");
0898   iDesc.add<std::string>("eid_output_name_energy", "output/regressed_energy");
0899   iDesc.add<std::string>("eid_output_name_id", "output/id_probabilities");
0900   iDesc.add<double>("eid_min_cluster_energy", 1.);
0901   iDesc.add<int>("eid_n_layers", 50);
0902   iDesc.add<int>("eid_n_clusters", 10);
0903   iDesc.add<bool>("computeLocalTime", false);
0904 }
0905 
0906 template class ticl::PatternRecognitionbyCLUE3D<TICLLayerTiles>;
0907 template class ticl::PatternRecognitionbyCLUE3D<TICLLayerTilesHFNose>;
[ 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!