Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2022-04-18 02:06:40

0001 #include "RecoLocalCalo/HGCalRecProducers/plugins/HGCalCLUEAlgo.h"
0002 
0003 // Geometry
0004 #include "DataFormats/HcalDetId/interface/HcalSubdetector.h"
0005 #include "Geometry/CaloGeometry/interface/CaloCellGeometry.h"
0006 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
0007 #include "Geometry/Records/interface/IdealGeometryRecord.h"
0008 
0009 #include "RecoEcal/EgammaCoreTools/interface/PositionCalc.h"
0010 //
0011 #include "DataFormats/CaloRecHit/interface/CaloID.h"
0012 #include "oneapi/tbb/task_arena.h"
0013 #include "oneapi/tbb.h"
0014 #include <limits>
0015 
0016 using namespace hgcal_clustering;
0017 
0018 template <typename T>
0019 void HGCalCLUEAlgoT<T>::getEventSetupPerAlgorithm(const edm::EventSetup& es) {
0020   cells_.clear();
0021   numberOfClustersPerLayer_.clear();
0022   cells_.resize(2 * (maxlayer_ + 1));
0023   numberOfClustersPerLayer_.resize(2 * (maxlayer_ + 1), 0);
0024 }
0025 
0026 template <typename T>
0027 void HGCalCLUEAlgoT<T>::populate(const HGCRecHitCollection& hits) {
0028   // loop over all hits and create the Hexel structure, skip energies below ecut
0029 
0030   if (dependSensor_) {
0031     // for each layer and wafer calculate the thresholds (sigmaNoise and energy)
0032     // once
0033     computeThreshold();
0034   }
0035 
0036   for (unsigned int i = 0; i < hits.size(); ++i) {
0037     const HGCRecHit& hgrh = hits[i];
0038     DetId detid = hgrh.detid();
0039     unsigned int layerOnSide = (rhtools_.getLayerWithOffset(detid) - 1);
0040 
0041     // set sigmaNoise default value 1 to use kappa value directly in case of
0042     // sensor-independent thresholds
0043     float sigmaNoise = 1.f;
0044     if (dependSensor_) {
0045       int thickness_index = rhtools_.getSiThickIndex(detid);
0046       if (thickness_index == -1)
0047         thickness_index = maxNumberOfThickIndices_;
0048 
0049       double storedThreshold = thresholds_[layerOnSide][thickness_index];
0050       if (detid.det() == DetId::HGCalHSi || detid.subdetId() == HGCHEF) {
0051         storedThreshold = thresholds_[layerOnSide][thickness_index + deltasi_index_regemfac_];
0052       }
0053       sigmaNoise = v_sigmaNoise_[layerOnSide][thickness_index];
0054 
0055       if (hgrh.energy() < storedThreshold)
0056         continue;  // this sets the ZS threshold at ecut times the sigma noise
0057                    // for the sensor
0058     }
0059     if (!dependSensor_ && hgrh.energy() < ecut_)
0060       continue;
0061     const GlobalPoint position(rhtools_.getPosition(detid));
0062     int offset = ((rhtools_.zside(detid) + 1) >> 1) * maxlayer_;
0063     int layer = layerOnSide + offset;
0064 
0065     cells_[layer].detid.emplace_back(detid);
0066     cells_[layer].x.emplace_back(position.x());
0067     cells_[layer].y.emplace_back(position.y());
0068     if (!rhtools_.isOnlySilicon(layer)) {
0069       cells_[layer].isSi.emplace_back(rhtools_.isSilicon(detid));
0070       cells_[layer].eta.emplace_back(position.eta());
0071       cells_[layer].phi.emplace_back(position.phi());
0072     }  // else, isSilicon == true and eta phi values will not be used
0073     cells_[layer].weight.emplace_back(hgrh.energy());
0074     cells_[layer].sigmaNoise.emplace_back(sigmaNoise);
0075   }
0076 }
0077 
0078 template <typename T>
0079 void HGCalCLUEAlgoT<T>::prepareDataStructures(unsigned int l) {
0080   auto cellsSize = cells_[l].detid.size();
0081   cells_[l].rho.resize(cellsSize, 0.f);
0082   cells_[l].delta.resize(cellsSize, 9999999);
0083   cells_[l].nearestHigher.resize(cellsSize, -1);
0084   cells_[l].clusterIndex.resize(cellsSize, -1);
0085   cells_[l].followers.resize(cellsSize);
0086   cells_[l].isSeed.resize(cellsSize, false);
0087   if (rhtools_.isOnlySilicon(l)) {
0088     cells_[l].isSi.resize(cellsSize, true);
0089     cells_[l].eta.resize(cellsSize, 0.f);
0090     cells_[l].phi.resize(cellsSize, 0.f);
0091   }
0092 }
0093 
0094 // Create a vector of Hexels associated to one cluster from a collection of
0095 // HGCalRecHits - this can be used directly to make the final cluster list -
0096 // this method can be invoked multiple times for the same event with different
0097 // input (reset should be called between events)
0098 template <typename T>
0099 void HGCalCLUEAlgoT<T>::makeClusters() {
0100   // assign all hits in each layer to a cluster core
0101   tbb::this_task_arena::isolate([&] {
0102     tbb::parallel_for(size_t(0), size_t(2 * maxlayer_ + 2), [&](size_t i) {
0103       prepareDataStructures(i);
0104       T lt;
0105       lt.clear();
0106       lt.fill(cells_[i].x, cells_[i].y, cells_[i].eta, cells_[i].phi, cells_[i].isSi);
0107       float delta_c;  // maximum search distance (critical distance) for local
0108                       // density calculation
0109       if (i % maxlayer_ < lastLayerEE_)
0110         delta_c = vecDeltas_[0];
0111       else if (i % maxlayer_ < (firstLayerBH_ - 1))
0112         delta_c = vecDeltas_[1];
0113       else
0114         delta_c = vecDeltas_[2];
0115       float delta_r = vecDeltas_[3];
0116       LogDebug("HGCalCLUEAlgo") << "maxlayer: " << maxlayer_ << " lastLayerEE: " << lastLayerEE_
0117                                 << " firstLayerBH: " << firstLayerBH_ << "\n";
0118 
0119       calculateLocalDensity(lt, i, delta_c, delta_r);
0120       calculateDistanceToHigher(lt, i, delta_c, delta_r);
0121       numberOfClustersPerLayer_[i] = findAndAssignClusters(i, delta_c, delta_r);
0122     });
0123   });
0124   //Now that we have the density per point we can store it
0125   for (unsigned int i = 0; i < 2 * maxlayer_ + 2; ++i) {
0126     setDensity(i);
0127   }
0128 }
0129 
0130 template <typename T>
0131 std::vector<reco::BasicCluster> HGCalCLUEAlgoT<T>::getClusters(bool) {
0132   std::vector<int> offsets(numberOfClustersPerLayer_.size(), 0);
0133 
0134   int maxClustersOnLayer = numberOfClustersPerLayer_[0];
0135 
0136   for (unsigned layerId = 1; layerId < offsets.size(); ++layerId) {
0137     offsets[layerId] = offsets[layerId - 1] + numberOfClustersPerLayer_[layerId - 1];
0138 
0139     maxClustersOnLayer = std::max(maxClustersOnLayer, numberOfClustersPerLayer_[layerId]);
0140   }
0141 
0142   auto totalNumberOfClusters = offsets.back() + numberOfClustersPerLayer_.back();
0143   clusters_v_.resize(totalNumberOfClusters);
0144   std::vector<std::vector<int>> cellsIdInCluster;
0145   cellsIdInCluster.reserve(maxClustersOnLayer);
0146 
0147   for (unsigned int layerId = 0; layerId < 2 * maxlayer_ + 2; ++layerId) {
0148     cellsIdInCluster.resize(numberOfClustersPerLayer_[layerId]);
0149     auto& cellsOnLayer = cells_[layerId];
0150     unsigned int numberOfCells = cellsOnLayer.detid.size();
0151     auto firstClusterIdx = offsets[layerId];
0152 
0153     for (unsigned int i = 0; i < numberOfCells; ++i) {
0154       auto clusterIndex = cellsOnLayer.clusterIndex[i];
0155       if (clusterIndex != -1)
0156         cellsIdInCluster[clusterIndex].push_back(i);
0157     }
0158 
0159     std::vector<std::pair<DetId, float>> thisCluster;
0160 
0161     for (auto& cl : cellsIdInCluster) {
0162       auto position = calculatePosition(cl, layerId);
0163       float energy = 0.f;
0164       int seedDetId = -1;
0165 
0166       for (auto cellIdx : cl) {
0167         energy += cellsOnLayer.weight[cellIdx];
0168         thisCluster.emplace_back(cellsOnLayer.detid[cellIdx], 1.f);
0169         if (cellsOnLayer.isSeed[cellIdx]) {
0170           seedDetId = cellsOnLayer.detid[cellIdx];
0171         }
0172       }
0173       auto globalClusterIndex = cellsOnLayer.clusterIndex[cl[0]] + firstClusterIdx;
0174 
0175       clusters_v_[globalClusterIndex] =
0176           reco::BasicCluster(energy, position, reco::CaloID::DET_HGCAL_ENDCAP, thisCluster, algoId_);
0177       clusters_v_[globalClusterIndex].setSeed(seedDetId);
0178       thisCluster.clear();
0179     }
0180 
0181     cellsIdInCluster.clear();
0182   }
0183   return clusters_v_;
0184 }
0185 
0186 template <typename T>
0187 math::XYZPoint HGCalCLUEAlgoT<T>::calculatePosition(const std::vector<int>& v, const unsigned int layerId) const {
0188   float total_weight = 0.f;
0189   float x = 0.f;
0190   float y = 0.f;
0191 
0192   unsigned int maxEnergyIndex = 0;
0193   float maxEnergyValue = 0.f;
0194 
0195   auto& cellsOnLayer = cells_[layerId];
0196 
0197   // loop over hits in cluster candidate
0198   // determining the maximum energy hit
0199   for (auto i : v) {
0200     total_weight += cellsOnLayer.weight[i];
0201     if (cellsOnLayer.weight[i] > maxEnergyValue) {
0202       maxEnergyValue = cellsOnLayer.weight[i];
0203       maxEnergyIndex = i;
0204     }
0205   }
0206 
0207   // Si cell or Scintillator. Used to set approach and parameters
0208   auto thick = rhtools_.getSiThickIndex(cellsOnLayer.detid[maxEnergyIndex]);
0209   bool isSiliconCell = (thick != -1);
0210 
0211   // TODO: this is recomputing everything twice and overwriting the position with log weighting position
0212   if (isSiliconCell) {
0213     float total_weight_log = 0.f;
0214     float x_log = 0.f;
0215     float y_log = 0.f;
0216     for (auto i : v) {
0217       //for silicon only just use 1+6 cells = 1.3cm for all thicknesses
0218       if (distance2(i, maxEnergyIndex, layerId, false) > positionDeltaRho2_)
0219         continue;
0220       float rhEnergy = cellsOnLayer.weight[i];
0221       float Wi = std::max(thresholdW0_[thick] + std::log(rhEnergy / total_weight), 0.);
0222       x_log += cellsOnLayer.x[i] * Wi;
0223       y_log += cellsOnLayer.y[i] * Wi;
0224       total_weight_log += Wi;
0225     }
0226 
0227     total_weight = total_weight_log;
0228     x = x_log;
0229     y = y_log;
0230   } else {
0231     for (auto i : v) {
0232       float rhEnergy = cellsOnLayer.weight[i];
0233 
0234       x += cellsOnLayer.x[i] * rhEnergy;
0235       y += cellsOnLayer.y[i] * rhEnergy;
0236     }
0237   }
0238   if (total_weight != 0.) {
0239     float inv_tot_weight = 1.f / total_weight;
0240     return math::XYZPoint(
0241         x * inv_tot_weight, y * inv_tot_weight, rhtools_.getPosition(cellsOnLayer.detid[maxEnergyIndex]).z());
0242   } else
0243     return math::XYZPoint(0.f, 0.f, 0.f);
0244 }
0245 
0246 template <typename T>
0247 void HGCalCLUEAlgoT<T>::calculateLocalDensity(const T& lt, const unsigned int layerId, float delta_c, float delta_r) {
0248   auto& cellsOnLayer = cells_[layerId];
0249   unsigned int numberOfCells = cellsOnLayer.detid.size();
0250   bool isOnlySi(false);
0251   if (rhtools_.isOnlySilicon(layerId))
0252     isOnlySi = true;
0253 
0254   for (unsigned int i = 0; i < numberOfCells; i++) {
0255     bool isSi = isOnlySi || cellsOnLayer.isSi[i];
0256     if (isSi) {
0257       float delta = delta_c;
0258       std::array<int, 4> search_box = lt.searchBox(
0259           cellsOnLayer.x[i] - delta, cellsOnLayer.x[i] + delta, cellsOnLayer.y[i] - delta, cellsOnLayer.y[i] + delta);
0260 
0261       for (int xBin = search_box[0]; xBin < search_box[1] + 1; ++xBin) {
0262         for (int yBin = search_box[2]; yBin < search_box[3] + 1; ++yBin) {
0263           int binId = lt.getGlobalBinByBin(xBin, yBin);
0264           size_t binSize = lt[binId].size();
0265 
0266           for (unsigned int j = 0; j < binSize; j++) {
0267             unsigned int otherId = lt[binId][j];
0268             bool otherSi = isOnlySi || cellsOnLayer.isSi[otherId];
0269             if (otherSi) {  //silicon cells cannot talk to scintillator cells
0270               if (distance(i, otherId, layerId, false) < delta) {
0271                 cellsOnLayer.rho[i] += (i == otherId ? 1.f : 0.5f) * cellsOnLayer.weight[otherId];
0272               }
0273             }
0274           }
0275         }
0276       }
0277     } else {
0278       float delta = delta_r;
0279       std::array<int, 4> search_box = lt.searchBoxEtaPhi(cellsOnLayer.eta[i] - delta,
0280                                                          cellsOnLayer.eta[i] + delta,
0281                                                          cellsOnLayer.phi[i] - delta,
0282                                                          cellsOnLayer.phi[i] + delta);
0283       cellsOnLayer.rho[i] += cellsOnLayer.weight[i];
0284       float northeast(0), northwest(0), southeast(0), southwest(0), all(0);
0285 
0286       for (int etaBin = search_box[0]; etaBin < search_box[1] + 1; ++etaBin) {
0287         for (int phiBin = search_box[2]; phiBin < search_box[3] + 1; ++phiBin) {
0288           int phi = (phiBin % T::type::nRowsPhi);
0289           int binId = lt.getGlobalBinByBinEtaPhi(etaBin, phi);
0290           size_t binSize = lt[binId].size();
0291 
0292           for (unsigned int j = 0; j < binSize; j++) {
0293             unsigned int otherId = lt[binId][j];
0294             if (!cellsOnLayer.isSi[otherId]) {  //scintillator cells cannot talk to silicon cells
0295               if (distance(i, otherId, layerId, true) < delta) {
0296                 int iPhi = HGCScintillatorDetId(cellsOnLayer.detid[i]).iphi();
0297                 int otherIPhi = HGCScintillatorDetId(cellsOnLayer.detid[otherId]).iphi();
0298                 int iEta = HGCScintillatorDetId(cellsOnLayer.detid[i]).ieta();
0299                 int otherIEta = HGCScintillatorDetId(cellsOnLayer.detid[otherId]).ieta();
0300                 int dIPhi = otherIPhi - iPhi;
0301                 dIPhi += abs(dIPhi) < 2 ? 0
0302                          : dIPhi < 0    ? scintMaxIphi_
0303                                         : -scintMaxIphi_;  // cells with iPhi=288 and iPhi=1 should be neiboring cells
0304                 int dIEta = otherIEta - iEta;
0305                 LogDebug("HGCalCLUEAlgo") << "  Debugging calculateLocalDensity for Scintillator: \n"
0306                                           << "    cell: " << otherId << " energy: " << cellsOnLayer.weight[otherId]
0307                                           << " otherIPhi: " << otherIPhi << " iPhi: " << iPhi
0308                                           << " otherIEta: " << otherIEta << " iEta: " << iEta << "\n";
0309 
0310                 if (otherId != i) {
0311                   auto neighborCellContribution = 0.5f * cellsOnLayer.weight[otherId];
0312                   all += neighborCellContribution;
0313                   if (dIPhi >= 0 && dIEta >= 0)
0314                     northeast += neighborCellContribution;
0315                   if (dIPhi <= 0 && dIEta >= 0)
0316                     southeast += neighborCellContribution;
0317                   if (dIPhi >= 0 && dIEta <= 0)
0318                     northwest += neighborCellContribution;
0319                   if (dIPhi <= 0 && dIEta <= 0)
0320                     southwest += neighborCellContribution;
0321                 }
0322                 LogDebug("HGCalCLUEAlgo") << "  Debugging calculateLocalDensity for Scintillator: \n"
0323                                           << "    northeast: " << northeast << " southeast: " << southeast
0324                                           << " northwest: " << northwest << " southwest: " << southwest << "\n";
0325               }
0326             }
0327           }
0328         }
0329       }
0330       float neighborsval = (std::max(northeast, northwest) > std::max(southeast, southwest))
0331                                ? std::max(northeast, northwest)
0332                                : std::max(southeast, southwest);
0333       if (use2x2_)
0334         cellsOnLayer.rho[i] += neighborsval;
0335       else
0336         cellsOnLayer.rho[i] += all;
0337     }
0338     LogDebug("HGCalCLUEAlgo") << "Debugging calculateLocalDensity: \n"
0339                               << "  cell: " << i << " isSilicon: " << cellsOnLayer.isSi[i]
0340                               << " eta: " << cellsOnLayer.eta[i] << " phi: " << cellsOnLayer.phi[i]
0341                               << " energy: " << cellsOnLayer.weight[i] << " density: " << cellsOnLayer.rho[i] << "\n";
0342   }
0343 }
0344 
0345 template <typename T>
0346 void HGCalCLUEAlgoT<T>::calculateDistanceToHigher(const T& lt,
0347                                                   const unsigned int layerId,
0348                                                   float delta_c,
0349                                                   float delta_r) {
0350   auto& cellsOnLayer = cells_[layerId];
0351   unsigned int numberOfCells = cellsOnLayer.detid.size();
0352   bool isOnlySi(false);
0353   if (rhtools_.isOnlySilicon(layerId))
0354     isOnlySi = true;
0355 
0356   for (unsigned int i = 0; i < numberOfCells; i++) {
0357     bool isSi = isOnlySi || cellsOnLayer.isSi[i];
0358     // initialize delta and nearest higher for i
0359     float maxDelta = std::numeric_limits<float>::max();
0360     float i_delta = maxDelta;
0361     int i_nearestHigher = -1;
0362     if (isSi) {
0363       float delta = delta_c;
0364       // get search box for ith hit
0365       // guarantee to cover a range "outlierDeltaFactor_*delta_c"
0366       auto range = outlierDeltaFactor_ * delta;
0367       std::array<int, 4> search_box = lt.searchBox(
0368           cellsOnLayer.x[i] - range, cellsOnLayer.x[i] + range, cellsOnLayer.y[i] - range, cellsOnLayer.y[i] + range);
0369       // loop over all bins in the search box
0370       for (int xBin = search_box[0]; xBin < search_box[1] + 1; ++xBin) {
0371         for (int yBin = search_box[2]; yBin < search_box[3] + 1; ++yBin) {
0372           // get the id of this bin
0373           size_t binId = lt.getGlobalBinByBin(xBin, yBin);
0374           // get the size of this bin
0375           size_t binSize = lt[binId].size();
0376 
0377           // loop over all hits in this bin
0378           for (unsigned int j = 0; j < binSize; j++) {
0379             unsigned int otherId = lt[binId][j];
0380             bool otherSi = isOnlySi || cellsOnLayer.isSi[otherId];
0381             if (otherSi) {  //silicon cells cannot talk to scintillator cells
0382               float dist = distance(i, otherId, layerId, false);
0383               bool foundHigher = (cellsOnLayer.rho[otherId] > cellsOnLayer.rho[i]) ||
0384                                  (cellsOnLayer.rho[otherId] == cellsOnLayer.rho[i] &&
0385                                   cellsOnLayer.detid[otherId] > cellsOnLayer.detid[i]);
0386               // if dist == i_delta, then last comer being the nearest higher
0387               if (foundHigher && dist <= i_delta) {
0388                 // update i_delta
0389                 i_delta = dist;
0390                 // update i_nearestHigher
0391                 i_nearestHigher = otherId;
0392               }
0393             }
0394           }
0395         }
0396       }
0397 
0398       bool foundNearestHigherInSearchBox = (i_delta != maxDelta);
0399       if (foundNearestHigherInSearchBox) {
0400         cellsOnLayer.delta[i] = i_delta;
0401         cellsOnLayer.nearestHigher[i] = i_nearestHigher;
0402       } else {
0403         // otherwise delta is guaranteed to be larger outlierDeltaFactor_*delta_c
0404         // we can safely maximize delta to be maxDelta
0405         cellsOnLayer.delta[i] = maxDelta;
0406         cellsOnLayer.nearestHigher[i] = -1;
0407       }
0408     } else {
0409       //similar to silicon
0410       float delta = delta_r;
0411       auto range = outlierDeltaFactor_ * delta;
0412       std::array<int, 4> search_box = lt.searchBoxEtaPhi(cellsOnLayer.eta[i] - range,
0413                                                          cellsOnLayer.eta[i] + range,
0414                                                          cellsOnLayer.phi[i] - range,
0415                                                          cellsOnLayer.phi[i] + range);
0416       // loop over all bins in the search box
0417       for (int xBin = search_box[0]; xBin < search_box[1] + 1; ++xBin) {
0418         for (int yBin = search_box[2]; yBin < search_box[3] + 1; ++yBin) {
0419           // get the id of this bin
0420           int phi = (yBin % T::type::nRowsPhi);
0421           size_t binId = lt.getGlobalBinByBinEtaPhi(xBin, phi);
0422           // get the size of this bin
0423           size_t binSize = lt[binId].size();
0424 
0425           // loop over all hits in this bin
0426           for (unsigned int j = 0; j < binSize; j++) {
0427             unsigned int otherId = lt[binId][j];
0428             if (!cellsOnLayer.isSi[otherId]) {  //scintillator cells cannot talk to silicon cells
0429               float dist = distance(i, otherId, layerId, true);
0430               bool foundHigher = (cellsOnLayer.rho[otherId] > cellsOnLayer.rho[i]) ||
0431                                  (cellsOnLayer.rho[otherId] == cellsOnLayer.rho[i] &&
0432                                   cellsOnLayer.detid[otherId] > cellsOnLayer.detid[i]);
0433               // if dist == i_delta, then last comer being the nearest higher
0434               if (foundHigher && dist <= i_delta) {
0435                 // update i_delta
0436                 i_delta = dist;
0437                 // update i_nearestHigher
0438                 i_nearestHigher = otherId;
0439               }
0440             }
0441           }
0442         }
0443       }
0444 
0445       bool foundNearestHigherInSearchBox = (i_delta != maxDelta);
0446       if (foundNearestHigherInSearchBox) {
0447         cellsOnLayer.delta[i] = i_delta;
0448         cellsOnLayer.nearestHigher[i] = i_nearestHigher;
0449       } else {
0450         // otherwise delta is guaranteed to be larger outlierDeltaFactor_*delta_c
0451         // we can safely maximize delta to be maxDelta
0452         cellsOnLayer.delta[i] = maxDelta;
0453         cellsOnLayer.nearestHigher[i] = -1;
0454       }
0455     }
0456     LogDebug("HGCalCLUEAlgo") << "Debugging calculateDistanceToHigher: \n"
0457                               << "  cell: " << i << " isSilicon: " << cellsOnLayer.isSi[i]
0458                               << " eta: " << cellsOnLayer.eta[i] << " phi: " << cellsOnLayer.phi[i]
0459                               << " energy: " << cellsOnLayer.weight[i] << " density: " << cellsOnLayer.rho[i]
0460                               << " nearest higher: " << cellsOnLayer.nearestHigher[i]
0461                               << " distance: " << cellsOnLayer.delta[i] << "\n";
0462   }
0463 }
0464 
0465 template <typename T>
0466 int HGCalCLUEAlgoT<T>::findAndAssignClusters(const unsigned int layerId, float delta_c, float delta_r) {
0467   // this is called once per layer and endcap...
0468   // so when filling the cluster temporary vector of Hexels we resize each time
0469   // by the number  of clusters found. This is always equal to the number of
0470   // cluster centers...
0471   unsigned int nClustersOnLayer = 0;
0472   auto& cellsOnLayer = cells_[layerId];
0473   unsigned int numberOfCells = cellsOnLayer.detid.size();
0474   std::vector<int> localStack;
0475   // find cluster seeds and outlier
0476   for (unsigned int i = 0; i < numberOfCells; i++) {
0477     float rho_c = kappa_ * cellsOnLayer.sigmaNoise[i];
0478     bool isSi = rhtools_.isOnlySilicon(layerId) || cellsOnLayer.isSi[i];
0479     float delta = isSi ? delta_c : delta_r;
0480 
0481     // initialize clusterIndex
0482     cellsOnLayer.clusterIndex[i] = -1;
0483     bool isSeed = (cellsOnLayer.delta[i] > delta) && (cellsOnLayer.rho[i] >= rho_c);
0484     bool isOutlier = (cellsOnLayer.delta[i] > outlierDeltaFactor_ * delta) && (cellsOnLayer.rho[i] < rho_c);
0485     if (isSeed) {
0486       cellsOnLayer.clusterIndex[i] = nClustersOnLayer;
0487       cellsOnLayer.isSeed[i] = true;
0488       nClustersOnLayer++;
0489       localStack.push_back(i);
0490 
0491     } else if (!isOutlier) {
0492       cellsOnLayer.followers[cellsOnLayer.nearestHigher[i]].push_back(i);
0493     }
0494   }
0495 
0496   // need to pass clusterIndex to their followers
0497   while (!localStack.empty()) {
0498     int endStack = localStack.back();
0499     auto& thisSeed = cellsOnLayer.followers[endStack];
0500     localStack.pop_back();
0501 
0502     // loop over followers
0503     for (int j : thisSeed) {
0504       // pass id to a follower
0505       cellsOnLayer.clusterIndex[j] = cellsOnLayer.clusterIndex[endStack];
0506       // push this follower to localStack
0507       localStack.push_back(j);
0508     }
0509   }
0510   return nClustersOnLayer;
0511 }
0512 
0513 template <typename T>
0514 void HGCalCLUEAlgoT<T>::computeThreshold() {
0515   // To support the TDR geometry and also the post-TDR one (v9 onwards), we
0516   // need to change the logic of the vectors containing signal to noise and
0517   // thresholds. The first 3 indices will keep on addressing the different
0518   // thicknesses of the Silicon detectors in CE_E , the next 3 indices will address
0519   // the thicknesses of the Silicon detectors in CE_H, while the last one, number 6 (the
0520   // seventh) will address the Scintillators. This change will support both
0521   // geometries at the same time.
0522 
0523   if (initialized_)
0524     return;  // only need to calculate thresholds once
0525 
0526   initialized_ = true;
0527 
0528   std::vector<double> dummy;
0529 
0530   dummy.resize(maxNumberOfThickIndices_ + !isNose_, 0);  // +1 to accomodate for the Scintillators
0531   thresholds_.resize(maxlayer_, dummy);
0532   v_sigmaNoise_.resize(maxlayer_, dummy);
0533 
0534   for (unsigned ilayer = 1; ilayer <= maxlayer_; ++ilayer) {
0535     for (unsigned ithick = 0; ithick < maxNumberOfThickIndices_; ++ithick) {
0536       float sigmaNoise = 0.001f * fcPerEle_ * nonAgedNoises_[ithick] * dEdXweights_[ilayer] /
0537                          (fcPerMip_[ithick] * thicknessCorrection_[ithick]);
0538       thresholds_[ilayer - 1][ithick] = sigmaNoise * ecut_;
0539       v_sigmaNoise_[ilayer - 1][ithick] = sigmaNoise;
0540       LogDebug("HGCalCLUEAlgo") << "ilayer: " << ilayer << " nonAgedNoises: " << nonAgedNoises_[ithick]
0541                                 << " fcPerEle: " << fcPerEle_ << " fcPerMip: " << fcPerMip_[ithick]
0542                                 << " noiseMip: " << fcPerEle_ * nonAgedNoises_[ithick] / fcPerMip_[ithick]
0543                                 << " sigmaNoise: " << sigmaNoise << "\n";
0544     }
0545 
0546     if (!isNose_) {
0547       float scintillators_sigmaNoise = 0.001f * noiseMip_ * dEdXweights_[ilayer] / sciThicknessCorrection_;
0548       thresholds_[ilayer - 1][maxNumberOfThickIndices_] = ecut_ * scintillators_sigmaNoise;
0549       v_sigmaNoise_[ilayer - 1][maxNumberOfThickIndices_] = scintillators_sigmaNoise;
0550       LogDebug("HGCalCLUEAlgo") << "ilayer: " << ilayer << " noiseMip: " << noiseMip_
0551                                 << " scintillators_sigmaNoise: " << scintillators_sigmaNoise << "\n";
0552     }
0553   }
0554 }
0555 
0556 template <typename T>
0557 void HGCalCLUEAlgoT<T>::setDensity(const unsigned int layerId) {
0558   auto& cellsOnLayer = cells_[layerId];
0559   unsigned int numberOfCells = cellsOnLayer.detid.size();
0560   for (unsigned int i = 0; i < numberOfCells; ++i)
0561     density_[cellsOnLayer.detid[i]] = cellsOnLayer.rho[i];
0562 }
0563 
0564 template <typename T>
0565 Density HGCalCLUEAlgoT<T>::getDensity() {
0566   return density_;
0567 }
0568 
0569 // explicit template instantiation
0570 template class HGCalCLUEAlgoT<HGCalLayerTiles>;
0571 template class HGCalCLUEAlgoT<HFNoseLayerTiles>;