Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2025-07-09 05:00:28 

0001 #include "RecoParticleFlow/PFClusterProducer/plugins/BarrelCLUEAlgo.h"
0002 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0003 
0004 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
0005 #include "DataFormats/CaloRecHit/interface/CaloID.h"
0006 #include "DataFormats/HcalDetId/interface/HcalDetId.h"
0007 
0008 #include "oneapi/tbb/task_arena.h"
0009 #include "oneapi/tbb.h"
0010 #include <limits>
0011 
0012 using namespace hgcal_clustering;
0013 
0014 template <typename T>
0015 void BarrelCLUEAlgoT<T>::setThresholds(edm::ESGetToken<EcalPFRecHitThresholds, EcalPFRecHitThresholdsRcd> EBThres,
0016                                        edm::ESGetToken<HcalPFCuts, HcalPFCutsRcd> HCALThres) {
0017   tok_ebThresholds_ = EBThres;
0018   tok_hcalThresholds_ = HCALThres;
0019 }
0020 
0021 template <typename T>
0022 void BarrelCLUEAlgoT<T>::getEventSetupPerAlgorithm(const edm::EventSetup& es) {
0023   cells_.clear();
0024   numberOfClustersPerLayer_.clear();
0025   ebThresholds_ = &es.getData(tok_ebThresholds_);
0026   hcalThresholds_ = &es.getData(tok_hcalThresholds_);
0027   maxlayer_ = maxLayerIndex_;  //for testing purposes
0028   cells_.resize(maxlayer_ + 1);
0029   numberOfClustersPerLayer_.resize(maxlayer_ + 1, 0);
0030 }
0031 
0032 template <typename T>
0033 void BarrelCLUEAlgoT<T>::populate(const reco::PFRecHitCollection& hits) {
0034   for (unsigned int i = 0; i < hits.size(); ++i) {
0035     const reco::PFRecHit& hit = hits[i];
0036     DetId detid(hit.detId());
0037     const GlobalPoint position = rhtools_.getPosition(detid);
0038     int layer = 0;
0039 
0040     if (detid.det() == DetId::Hcal) {
0041       HcalDetId hid(detid);
0042       layer = hid.depth();
0043       if (detid.subdetId() == HcalSubdetector::HcalOuter)
0044         layer += 1;
0045     }
0046 
0047     cells_[layer].detid.push_back(detid);
0048     cells_[layer].eta.push_back(position.eta());
0049     cells_[layer].phi.push_back(position.phi());
0050     cells_[layer].weight.push_back(hit.energy());
0051     cells_[layer].r.push_back(position.mag());
0052     float sigmaNoise = 0.f;
0053     if (detid.det() == DetId::Ecal) {
0054       sigmaNoise = (*ebThresholds_)[detid];
0055     } else {
0056       const HcalPFCut* item = (*hcalThresholds_).getValues(detid);
0057       sigmaNoise = item->seedThreshold();
0058     }
0059     cells_[layer].sigmaNoise.push_back(sigmaNoise);
0060   }
0061 }
0062 
0063 template <typename T>
0064 void BarrelCLUEAlgoT<T>::prepareDataStructures(unsigned int l) {
0065   auto cellsSize = cells_[l].detid.size();
0066   cells_[l].rho.resize(cellsSize, 0.f);
0067   cells_[l].delta.resize(cellsSize, 9999999);
0068   cells_[l].nearestHigher.resize(cellsSize, -1);
0069   cells_[l].clusterIndex.resize(cellsSize);
0070   cells_[l].followers.resize(cellsSize);
0071   cells_[l].isSeed.resize(cellsSize, false);
0072   cells_[l].eta.resize(cellsSize, 0.f);
0073   cells_[l].phi.resize(cellsSize, 0.f);
0074   cells_[l].r.resize(cellsSize, 0.f);
0075 }
0076 
0077 template <typename T>
0078 void BarrelCLUEAlgoT<T>::makeClusters() {
0079   tbb::this_task_arena::isolate([&] {
0080     tbb::parallel_for(size_t(0), size_t(maxlayer_ + 1), [&](size_t i) {
0081       prepareDataStructures(i);
0082       T lt;
0083       lt.clear();
0084       lt.fill(cells_[i].eta, cells_[i].phi);
0085 
0086       calculateLocalDensity(lt, i, deltac_);
0087       calculateDistanceToHigher(lt, i, deltac_);
0088       numberOfClustersPerLayer_[i] = findAndAssignClusters(i, deltac_);
0089 
0090       if (doSharing_)
0091         // Now running the sharing routine
0092         passSharedClusterIndex(lt, i, deltac_);
0093     });
0094   });
0095   for (unsigned int i = 0; i < maxlayer_ + 1; ++i) {
0096     setDensity(i);
0097   }
0098 }
0099 
0100 template <typename T>
0101 std::vector<reco::BasicCluster> BarrelCLUEAlgoT<T>::getClusters(bool) {
0102   std::vector<int> offsets(numberOfClustersPerLayer_.size(), 0);
0103 
0104   int maxClustersOnLayer = numberOfClustersPerLayer_[0];
0105 
0106   if constexpr (!std::is_same_v<T, EBLayerTiles>) {
0107     for (unsigned layerId = 1; layerId < offsets.size(); ++layerId) {
0108       offsets[layerId] = offsets[layerId - 1] + numberOfClustersPerLayer_[layerId - 1];
0109       maxClustersOnLayer = std::max(maxClustersOnLayer, numberOfClustersPerLayer_[layerId]);
0110     }
0111   }
0112 
0113   auto totalNumberOfClusters = offsets.back() + numberOfClustersPerLayer_.back();
0114   clusters_v_.resize(totalNumberOfClusters);
0115   // Store the cellId and energy for each cluster
0116   // We need to store the energy as a pair here because we perform the energy splitting
0117   std::vector<std::vector<std::pair<int, float>>> cellsIdInCluster;
0118   cellsIdInCluster.reserve(maxClustersOnLayer);
0119 
0120   for (unsigned int layerId = 0; layerId < maxlayer_ + 1; ++layerId) {
0121     cellsIdInCluster.resize(numberOfClustersPerLayer_[layerId]);
0122     auto& cellsOnLayer = cells_[layerId];
0123     unsigned int numberOfCells = cellsOnLayer.detid.size();
0124     auto firstClusterIdx = offsets[layerId];
0125 
0126     for (unsigned int i = 0; i < numberOfCells; ++i) {
0127       auto clusterIndex = cellsOnLayer.clusterIndex[i];
0128 
0129       if (clusterIndex.size() == 1) {
0130         cellsIdInCluster[clusterIndex[0]].push_back(std::make_pair(i, 1.f));
0131       } else if (clusterIndex.size() > 1) {
0132         std::vector<float> fractions(clusterIndex.size());
0133 
0134         for (unsigned int j = 0; j < clusterIndex.size(); j++) {
0135           const auto& seed = clusterIndex[j];
0136           const auto& seedCell = cellsOnLayer.seedToCellIndex[seed];
0137           const auto& seedEnergy = cellsOnLayer.weight[seedCell];
0138           // compute the distance
0139           float dist = distance(i, seedCell, layerId) / T::type::cellWidthEta;
0140           fractions[j] = seedEnergy * std::exp(-(std::pow(dist, 2)) / (2 * std::pow(T::type::showerSigma, 2)));
0141         }
0142         auto tot_norm_fractions = std::accumulate(std::begin(fractions), std::end(fractions), 0.);
0143 
0144         for (unsigned int j = 0; j < clusterIndex.size(); j++) {
0145           float norm_fraction = fractions[j] / tot_norm_fractions;
0146           if (norm_fraction < fractionCutoff_) {
0147             auto clusterIndex_it = std::find(clusterIndex.begin(), clusterIndex.end(), clusterIndex[j]);
0148             auto fraction_it = std::find(fractions.begin(), fractions.end(), fractions[j]);
0149             clusterIndex.erase(clusterIndex_it);
0150             fractions.erase(fraction_it);
0151           }
0152           if (norm_fraction >= fractionCutoff_) {
0153             cellsIdInCluster[clusterIndex[j]].push_back(std::make_pair(i, norm_fraction));
0154           }
0155         }
0156         auto tot_norm_cleaned_fractions = std::accumulate(fractions.begin(), fractions.end(), 0.);
0157         for (unsigned int j = 0; j < clusterIndex.size(); j++) {
0158           float norm_fraction = fractions[j] / tot_norm_cleaned_fractions;
0159           cellsIdInCluster[clusterIndex[j]].push_back(std::make_pair(i, norm_fraction));
0160         }
0161       }
0162     }
0163 
0164     std::vector<std::pair<DetId, float>> thisCluster;
0165     for (int clIndex = 0; clIndex < numberOfClustersPerLayer_[layerId]; clIndex++) {
0166       auto& cl = cellsIdInCluster[clIndex];
0167       auto position = calculatePosition(cl, layerId);
0168       float energy = 0.f;
0169       int seedDetId = -1;
0170 
0171       for (auto [cellIdx, fraction] : cl) {
0172         energy += cellsOnLayer.weight[cellIdx] * fraction;
0173         thisCluster.emplace_back(cellsOnLayer.detid[cellIdx], fraction);
0174         if (cellsOnLayer.isSeed[cellIdx]) {
0175           seedDetId = cellsOnLayer.detid[cellIdx];
0176         }
0177       }
0178       auto globalClusterIndex = clIndex + firstClusterIdx;
0179 
0180       if constexpr (std::is_same_v<T, EBLayerTiles>) {
0181         clusters_v_[globalClusterIndex] =
0182             reco::BasicCluster(energy, position, reco::CaloID::DET_ECAL_BARREL, thisCluster, algoId_);
0183       } else if constexpr (std::is_same_v<T, HBLayerTiles>) {
0184         clusters_v_[globalClusterIndex] =
0185             reco::BasicCluster(energy, position, reco::CaloID::DET_HCAL_BARREL, thisCluster, algoId_);
0186       } else {
0187         clusters_v_[globalClusterIndex] =
0188             reco::BasicCluster(energy, position, reco::CaloID::DET_HO, thisCluster, algoId_);
0189       }
0190       clusters_v_[globalClusterIndex].setSeed(seedDetId);
0191       thisCluster.clear();
0192     }
0193     cellsIdInCluster.clear();
0194   }
0195   return clusters_v_;
0196 }
0197 
0198 template <typename T>
0199 math::XYZPoint BarrelCLUEAlgoT<T>::calculatePosition(const std::vector<std::pair<int, float>>& v,
0200                                                      const unsigned int layerId) const {
0201   float total_weight = 0.f;
0202   float x = 0.f;
0203   float y = 0.f;
0204   float z = 0.f;
0205 
0206   auto& cellsOnLayer = cells_[layerId];
0207   for (const auto& [i, fraction] : v) {
0208     float rhEnergy = cellsOnLayer.weight[i] * fraction;
0209     total_weight += rhEnergy;
0210     float theta = 2 * std::atan(std::exp(-cellsOnLayer.eta[i]));
0211     const GlobalPoint gp(GlobalPoint::Polar(theta, cellsOnLayer.phi[i], cellsOnLayer.r[i]));
0212     x += gp.x() * rhEnergy;
0213     y += gp.y() * rhEnergy;
0214     z += gp.z() * rhEnergy;
0215   }
0216   if (total_weight != 0.f) {
0217     float inv_total_weight = 1.f / total_weight;
0218     return math::XYZPoint(x * inv_total_weight, y * inv_total_weight, z * inv_total_weight);
0219   } else {
0220     return math::XYZPoint(0.f, 0.f, 0.f);
0221   }
0222 }
0223 
0224 template <typename T>
0225 void BarrelCLUEAlgoT<T>::calculateLocalDensity(const T& lt, const unsigned int layerId, float deltac_) {
0226   auto& cellsOnLayer = cells_[layerId];
0227   unsigned int numberOfCells = cellsOnLayer.detid.size();
0228 
0229   for (unsigned int i = 0; i < numberOfCells; ++i) {
0230     float delta = deltac_;
0231     std::array<int, 4> search_box = lt.searchBox(cellsOnLayer.eta[i] - delta,
0232                                                  cellsOnLayer.eta[i] + delta,
0233                                                  cellsOnLayer.phi[i] - delta,
0234                                                  cellsOnLayer.phi[i] + delta);
0235     LogDebug("BarrelCLUEAlgo") << "searchBox for cell " << i << " cell eta: " << cellsOnLayer.eta[i]
0236                                << " cell phi: " << cellsOnLayer.phi[i] << " eta window: " << cellsOnLayer.eta[i] - delta
0237                                << " -> " << cellsOnLayer.eta[i] + delta
0238                                << " phi window: " << cellsOnLayer.phi[i] - delta << " -> "
0239                                << cellsOnLayer.phi[i] + delta << "\n";
0240     for (int etaBin = search_box[0]; etaBin < search_box[1]; ++etaBin) {
0241       for (int phiBin = search_box[2]; phiBin < search_box[3]; ++phiBin) {
0242         int phi = (phiBin % T::type::nRows);
0243         int binId = lt.getGlobalBinByBin(etaBin, phi);
0244         size_t binSize = lt[binId].size();
0245         LogDebug("BarrelCLUEAlgo") << "Bin size for cell " << i << ": " << binSize << " binId: " << binId
0246                                    << " etaBin: " << etaBin << " phi: " << phi << "\n";
0247         for (unsigned int j = 0; j < binSize; ++j) {
0248           unsigned int otherId = lt[binId][j];
0249           if (distance(i, otherId, layerId) < delta) {
0250             cellsOnLayer.rho[i] += (i == otherId ? 1.f : 0.5f) * cellsOnLayer.weight[otherId];
0251           }
0252         }
0253       }
0254     }
0255   }
0256 }
0257 
0258 template <typename T>
0259 void BarrelCLUEAlgoT<T>::calculateDistanceToHigher(const T& lt, const unsigned int layerId, float delta_c) {
0260   auto& cellsOnLayer = cells_[layerId];
0261   unsigned int numberOfCells = cellsOnLayer.detid.size();
0262 
0263   for (unsigned int i = 0; i < numberOfCells; ++i) {
0264     float maxDelta = std::numeric_limits<float>::max();
0265     float i_delta = maxDelta;
0266     float i_nearestHigher = -1;
0267 
0268     float range = delta_c;
0269     std::array<int, 4> search_box = lt.searchBox(cellsOnLayer.eta[i] - range,
0270                                                  cellsOnLayer.eta[i] + range,
0271                                                  cellsOnLayer.phi[i] - range,
0272                                                  cellsOnLayer.phi[i] + range);
0273 
0274     for (int etaBin = search_box[0]; etaBin < search_box[1]; ++etaBin) {
0275       for (int phiBin = search_box[2]; phiBin < search_box[3]; ++phiBin) {
0276         int phi = (phiBin % T::type::nRows);
0277         size_t binId = lt.getGlobalBinByBin(etaBin, phi);
0278         size_t binSize = lt[binId].size();
0279 
0280         for (unsigned int j = 0; j < binSize; ++j) {
0281           int otherId = lt[binId][j];
0282           float dist = distance(i, otherId, layerId);
0283           bool foundHigher =
0284               (cellsOnLayer.rho[otherId] > cellsOnLayer.rho[i]) ||
0285               (cellsOnLayer.rho[otherId] == cellsOnLayer.rho[i] && cellsOnLayer.detid[otherId] > cellsOnLayer.detid[i]);
0286           if (foundHigher && dist <= i_delta) {
0287             i_delta = dist;
0288             i_nearestHigher = otherId;
0289           }
0290         }
0291       }
0292     }
0293     bool foundNearestHigherInSearchBox = (i_delta != maxDelta);
0294     if (foundNearestHigherInSearchBox) {
0295       cellsOnLayer.delta[i] = i_delta;
0296       cellsOnLayer.nearestHigher[i] = i_nearestHigher;
0297     } else {
0298       cellsOnLayer.delta[i] = maxDelta;
0299       cellsOnLayer.nearestHigher[i] = -1;
0300     }
0301     LogDebug("BarrelCLUEAlgo") << "Debugging calculateDistanceToHigher: \n"
0302                                << " cell: " << i << " eta: " << cellsOnLayer.eta[i] << " phi: " << cellsOnLayer.phi[i]
0303                                << " energy: " << cellsOnLayer.weight[i] << " density: " << cellsOnLayer.rho[i]
0304                                << " nearest higher: " << cellsOnLayer.nearestHigher[i]
0305                                << " distance: " << cellsOnLayer.delta[i] << "\n";
0306   }
0307 }
0308 
0309 template <typename T>
0310 int BarrelCLUEAlgoT<T>::findAndAssignClusters(const unsigned int layerId, float delta_c) {
0311   unsigned int nClustersOnLayer = 0;
0312   auto& cellsOnLayer = cells_[layerId];
0313   unsigned int numberOfCells = cellsOnLayer.detid.size();
0314   std::vector<int> localStack;
0315   for (unsigned int i = 0; i < numberOfCells; ++i) {
0316     float rho_c = kappa_ * cellsOnLayer.sigmaNoise[i];
0317     float delta = delta_c;
0318 
0319     //std::vector<float> outDeltas = {7., 9., 10., 11.};
0320     //if constexpr (std::is_same_v<T, HBLayerTiles>)
0321     //  outlierDeltaFactor_ = outDeltas[layerId - 1];
0322     bool isSeed = (cellsOnLayer.delta[i] > delta) && (cellsOnLayer.rho[i] >= rho_c);
0323     bool isOutlier = (cellsOnLayer.delta[i] > outlierDeltaFactor_) && (cellsOnLayer.rho[i] < rho_c);
0324     if (isSeed) {
0325       cellsOnLayer.clusterIndex[i].push_back(nClustersOnLayer);
0326       cellsOnLayer.isSeed[i] = true;
0327       cellsOnLayer.seedToCellIndex.push_back(i);
0328       nClustersOnLayer++;
0329       localStack.push_back(i);
0330     } else if (!isOutlier) {
0331       cellsOnLayer.followers[cellsOnLayer.nearestHigher[i]].push_back(i);
0332     }
0333   }
0334 
0335   while (!localStack.empty()) {
0336     int endStack = localStack.back();
0337     auto& thisSeed = cellsOnLayer.followers[endStack];
0338     localStack.pop_back();
0339 
0340     for (int j : thisSeed) {
0341       cellsOnLayer.clusterIndex[j].push_back(cellsOnLayer.clusterIndex[endStack][0]);
0342       localStack.push_back(j);
0343     }
0344   }
0345   return nClustersOnLayer;
0346 }
0347 
0348 template <typename T>
0349 void BarrelCLUEAlgoT<T>::passSharedClusterIndex(const T& lt, const unsigned int layerId, float delta_c) {
0350   auto& cellsOnLayer = cells_[layerId];
0351   unsigned int numberOfCells = cellsOnLayer.detid.size();
0352   float delta = delta_c;
0353   for (unsigned int i = 0; i < numberOfCells; ++i) {
0354     // Do not run on outliers, but run also on seeds and followers
0355     if ((cellsOnLayer.clusterIndex[i].size() == 0))
0356       continue;
0357 
0358     std::array<int, 4> search_box = lt.searchBox(cellsOnLayer.eta[i] - delta,
0359                                                  cellsOnLayer.eta[i] + delta,
0360                                                  cellsOnLayer.phi[i] - delta,
0361                                                  cellsOnLayer.phi[i] + delta);
0362 
0363     for (int etaBin = search_box[0]; etaBin < search_box[1]; ++etaBin) {
0364       for (int phiBin = search_box[2]; phiBin < search_box[3]; ++phiBin) {
0365         int phi = (phiBin % T::type::nRows);
0366         size_t binId = lt.getGlobalBinByBin(etaBin, phi);
0367         size_t binSize = lt[binId].size();
0368 
0369         for (unsigned int j = 0; j < binSize; ++j) {
0370           unsigned int otherId = lt[binId][j];
0371           if (cellsOnLayer.clusterIndex[otherId].size() == 0)
0372             continue;
0373           int otherClusterIndex = cellsOnLayer.clusterIndex[otherId][0];
0374           if (std::find(std::begin(cellsOnLayer.clusterIndex[i]),
0375                         std::end(cellsOnLayer.clusterIndex[i]),
0376                         otherClusterIndex) == std::end(cellsOnLayer.clusterIndex[i]))
0377             cellsOnLayer.clusterIndex[i].push_back(otherClusterIndex);
0378         }
0379       }
0380     }
0381   }
0382 }
0383 
0384 template <typename T>
0385 void BarrelCLUEAlgoT<T>::setDensity(const unsigned int layerId) {
0386   auto& cellsOnLayer = cells_[layerId];
0387   unsigned int numberOfCells = cellsOnLayer.detid.size();
0388   for (unsigned int i = 0; i < numberOfCells; ++i) {
0389     density_[cellsOnLayer.detid[i]] = cellsOnLayer.rho[i];
0390   }
0391 }
0392 
0393 template <typename T>
0394 Density BarrelCLUEAlgoT<T>::getDensity() {
0395   return density_;
0396 }
0397 
0398 template class BarrelCLUEAlgoT<EBLayerTiles>;
0399 template class BarrelCLUEAlgoT<HBLayerTiles>;
[ 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!