Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2022-08-30 22:39:07 

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