File indexing completed on 2023-06-07 02:17:53
0001 #include "RecoLocalCalo/HGCalRecProducers/plugins/HGCalCLUEAlgo.h"
0002
0003
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 #include "DataFormats/DetId/interface/DetId.h"
0016
0017 using namespace hgcal_clustering;
0018
0019 template <typename T, typename STRATEGY>
0020 void HGCalCLUEAlgoT<T, STRATEGY>::getEventSetupPerAlgorithm(const edm::EventSetup& es) {
0021 cells_.clear();
0022 numberOfClustersPerLayer_.clear();
0023 cells_.resize(2 * (maxlayer_ + 1));
0024 numberOfClustersPerLayer_.resize(2 * (maxlayer_ + 1), 0);
0025 }
0026
0027 template <typename T, typename STRATEGY>
0028 void HGCalCLUEAlgoT<T, STRATEGY>::populate(const HGCRecHitCollection& hits) {
0029
0030 if (dependSensor_) {
0031
0032
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
0042
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;
0057
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 if constexpr (std::is_same_v<STRATEGY, HGCalScintillatorStrategy>) {
0067 cells_[layer].dim1.emplace_back(position.eta());
0068 cells_[layer].dim2.emplace_back(position.phi());
0069 }
0070 else {
0071 cells_[layer].dim1.emplace_back(position.x());
0072 cells_[layer].dim2.emplace_back(position.y());
0073 }
0074 cells_[layer].weight.emplace_back(hgrh.energy());
0075 cells_[layer].sigmaNoise.emplace_back(sigmaNoise);
0076 }
0077 }
0078
0079 template <typename T, typename STRATEGY>
0080 void HGCalCLUEAlgoT<T, STRATEGY>::prepareDataStructures(unsigned int l) {
0081 auto cellsSize = cells_[l].detid.size();
0082 cells_[l].rho.resize(cellsSize, 0.f);
0083 cells_[l].delta.resize(cellsSize, 9999999);
0084 cells_[l].nearestHigher.resize(cellsSize, -1);
0085 cells_[l].clusterIndex.resize(cellsSize, -1);
0086 cells_[l].followers.resize(cellsSize);
0087 cells_[l].isSeed.resize(cellsSize, false);
0088 }
0089
0090
0091
0092
0093
0094 template <typename T, typename STRATEGY>
0095 void HGCalCLUEAlgoT<T, STRATEGY>::makeClusters() {
0096
0097 tbb::this_task_arena::isolate([&] {
0098 tbb::parallel_for(size_t(0), size_t(2 * maxlayer_ + 2), [&](size_t i) {
0099 prepareDataStructures(i);
0100 T lt;
0101 lt.clear();
0102 lt.fill(cells_[i].dim1, cells_[i].dim2);
0103
0104 float delta;
0105 if constexpr (std::is_same_v<STRATEGY, HGCalSiliconStrategy>) {
0106
0107 float delta_c;
0108 if (i % maxlayer_ < lastLayerEE_)
0109 delta_c = vecDeltas_[0];
0110 else if (i % maxlayer_ < (firstLayerBH_ - 1))
0111 delta_c = vecDeltas_[1];
0112 else
0113 delta_c = vecDeltas_[2];
0114 delta = delta_c;
0115 } else {
0116 float delta_r = vecDeltas_[3];
0117 delta = delta_r;
0118 }
0119 LogDebug("HGCalCLUEAlgo") << "maxlayer: " << maxlayer_ << " lastLayerEE: " << lastLayerEE_
0120 << " firstLayerBH: " << firstLayerBH_ << "\n";
0121
0122 calculateLocalDensity(lt, i, delta);
0123 calculateDistanceToHigher(lt, i, delta);
0124 numberOfClustersPerLayer_[i] = findAndAssignClusters(i, delta);
0125 });
0126 });
0127 }
0128
0129 template <typename T, typename STRATEGY>
0130 std::vector<reco::BasicCluster> HGCalCLUEAlgoT<T, STRATEGY>::getClusters(bool) {
0131 std::vector<int> offsets(numberOfClustersPerLayer_.size(), 0);
0132
0133 int maxClustersOnLayer = numberOfClustersPerLayer_[0];
0134
0135 for (unsigned layerId = 1; layerId < offsets.size(); ++layerId) {
0136 offsets[layerId] = offsets[layerId - 1] + numberOfClustersPerLayer_[layerId - 1];
0137
0138 maxClustersOnLayer = std::max(maxClustersOnLayer, numberOfClustersPerLayer_[layerId]);
0139 }
0140
0141 auto totalNumberOfClusters = offsets.back() + numberOfClustersPerLayer_.back();
0142 clusters_v_.resize(totalNumberOfClusters);
0143 std::vector<std::vector<int>> cellsIdInCluster;
0144 cellsIdInCluster.reserve(maxClustersOnLayer);
0145
0146 for (unsigned int layerId = 0; layerId < 2 * maxlayer_ + 2; ++layerId) {
0147 cellsIdInCluster.resize(numberOfClustersPerLayer_[layerId]);
0148 auto& cellsOnLayer = cells_[layerId];
0149 unsigned int numberOfCells = cellsOnLayer.detid.size();
0150 auto firstClusterIdx = offsets[layerId];
0151
0152 for (unsigned int i = 0; i < numberOfCells; ++i) {
0153 auto clusterIndex = cellsOnLayer.clusterIndex[i];
0154 if (clusterIndex != -1)
0155 cellsIdInCluster[clusterIndex].push_back(i);
0156 }
0157
0158 std::vector<std::pair<DetId, float>> thisCluster;
0159
0160 for (auto& cl : cellsIdInCluster) {
0161 math::XYZPoint position = math::XYZPoint(0.f, 0.f, 0.f);
0162 float energy = 0.f;
0163 int seedDetId = -1;
0164
0165 for (auto cellIdx : cl) {
0166 energy += cellsOnLayer.weight[cellIdx];
0167 thisCluster.emplace_back(cellsOnLayer.detid[cellIdx], 1.f);
0168 if (cellsOnLayer.isSeed[cellIdx]) {
0169 seedDetId = cellsOnLayer.detid[cellIdx];
0170 }
0171 }
0172 auto globalClusterIndex = cellsOnLayer.clusterIndex[cl[0]] + firstClusterIdx;
0173
0174 clusters_v_[globalClusterIndex] =
0175 reco::BasicCluster(energy, position, reco::CaloID::DET_HGCAL_ENDCAP, thisCluster, algoId_);
0176 clusters_v_[globalClusterIndex].setSeed(seedDetId);
0177 thisCluster.clear();
0178 }
0179
0180 cellsIdInCluster.clear();
0181 }
0182 return clusters_v_;
0183 }
0184 template <typename T, typename STRATEGY>
0185 void HGCalCLUEAlgoT<T, STRATEGY>::calculateLocalDensity(const T& lt,
0186 const unsigned int layerId,
0187 float delta,
0188 HGCalSiliconStrategy strategy) {
0189 auto& cellsOnLayer = cells_[layerId];
0190 unsigned int numberOfCells = cellsOnLayer.detid.size();
0191 for (unsigned int i = 0; i < numberOfCells; i++) {
0192 std::array<int, 4> search_box = lt.searchBox(cellsOnLayer.dim1[i] - delta,
0193 cellsOnLayer.dim1[i] + delta,
0194 cellsOnLayer.dim2[i] - delta,
0195 cellsOnLayer.dim2[i] + delta);
0196
0197 for (int xBin = search_box[0]; xBin < search_box[1] + 1; ++xBin) {
0198 for (int yBin = search_box[2]; yBin < search_box[3] + 1; ++yBin) {
0199 int binId = lt.getGlobalBinByBin(xBin, yBin);
0200 size_t binSize = lt[binId].size();
0201
0202 for (unsigned int j = 0; j < binSize; j++) {
0203 unsigned int otherId = lt[binId][j];
0204 if (distance(lt, i, otherId, layerId) < delta) {
0205 cellsOnLayer.rho[i] += (i == otherId ? 1.f : 0.5f) * cellsOnLayer.weight[otherId];
0206 }
0207 }
0208 }
0209 }
0210 LogDebug("HGCalCLUEAlgo") << "Debugging calculateLocalDensity: \n"
0211 << " cell: " << i << " eta: " << cellsOnLayer.dim1[i] << " phi: " << cellsOnLayer.dim2[i]
0212 << " energy: " << cellsOnLayer.weight[i] << " density: " << cellsOnLayer.rho[i] << "\n";
0213 }
0214 }
0215 template <typename T, typename STRATEGY>
0216 void HGCalCLUEAlgoT<T, STRATEGY>::calculateLocalDensity(const T& lt,
0217 const unsigned int layerId,
0218 float delta,
0219 HGCalScintillatorStrategy strategy) {
0220 auto& cellsOnLayer = cells_[layerId];
0221 unsigned int numberOfCells = cellsOnLayer.detid.size();
0222 for (unsigned int i = 0; i < numberOfCells; i++) {
0223 std::array<int, 4> search_box = lt.searchBox(cellsOnLayer.dim1[i] - delta,
0224 cellsOnLayer.dim1[i] + delta,
0225 cellsOnLayer.dim2[i] - delta,
0226 cellsOnLayer.dim2[i] + delta);
0227 cellsOnLayer.rho[i] += cellsOnLayer.weight[i];
0228 float northeast(0), northwest(0), southeast(0), southwest(0), all(0);
0229 for (int etaBin = search_box[0]; etaBin < search_box[1] + 1; ++etaBin) {
0230 for (int phiBin = search_box[2]; phiBin < search_box[3] + 1; ++phiBin) {
0231 int phi = (phiBin % T::type::nRows);
0232 int binId = lt.getGlobalBinByBin(etaBin, phi);
0233 size_t binSize = lt[binId].size();
0234 for (unsigned int j = 0; j < binSize; j++) {
0235 unsigned int otherId = lt[binId][j];
0236 if (distance(lt, i, otherId, layerId) < delta) {
0237 int iPhi = HGCScintillatorDetId(cellsOnLayer.detid[i]).iphi();
0238 int otherIPhi = HGCScintillatorDetId(cellsOnLayer.detid[otherId]).iphi();
0239 int iEta = HGCScintillatorDetId(cellsOnLayer.detid[i]).ieta();
0240 int otherIEta = HGCScintillatorDetId(cellsOnLayer.detid[otherId]).ieta();
0241 int dIPhi = otherIPhi - iPhi;
0242 dIPhi += abs(dIPhi) < 2 ? 0
0243 : dIPhi < 0 ? scintMaxIphi_
0244 : -scintMaxIphi_;
0245 int dIEta = otherIEta - iEta;
0246 LogDebug("HGCalCLUEAlgo") << " Debugging calculateLocalDensity for Scintillator: \n"
0247 << " cell: " << otherId << " energy: " << cellsOnLayer.weight[otherId]
0248 << " otherIPhi: " << otherIPhi << " iPhi: " << iPhi << " otherIEta: " << otherIEta
0249 << " iEta: " << iEta << "\n";
0250
0251 if (otherId != i) {
0252 auto neighborCellContribution = 0.5f * cellsOnLayer.weight[otherId];
0253 all += neighborCellContribution;
0254 if (dIPhi >= 0 && dIEta >= 0)
0255 northeast += neighborCellContribution;
0256 if (dIPhi <= 0 && dIEta >= 0)
0257 southeast += neighborCellContribution;
0258 if (dIPhi >= 0 && dIEta <= 0)
0259 northwest += neighborCellContribution;
0260 if (dIPhi <= 0 && dIEta <= 0)
0261 southwest += neighborCellContribution;
0262 }
0263 LogDebug("HGCalCLUEAlgo") << " Debugging calculateLocalDensity for Scintillator: \n"
0264 << " northeast: " << northeast << " southeast: " << southeast
0265 << " northwest: " << northwest << " southwest: " << southwest << "\n";
0266 }
0267 }
0268 }
0269 }
0270 float neighborsval = (std::max(northeast, northwest) > std::max(southeast, southwest))
0271 ? std::max(northeast, northwest)
0272 : std::max(southeast, southwest);
0273 if (use2x2_)
0274 cellsOnLayer.rho[i] += neighborsval;
0275 else
0276 cellsOnLayer.rho[i] += all;
0277 LogDebug("HGCalCLUEAlgo") << "Debugging calculateLocalDensity: \n"
0278 << " cell: " << i << " eta: " << cellsOnLayer.dim1[i] << " phi: " << cellsOnLayer.dim2[i]
0279 << " energy: " << cellsOnLayer.weight[i] << " density: " << cellsOnLayer.rho[i] << "\n";
0280 }
0281 }
0282 template <typename T, typename STRATEGY>
0283 void HGCalCLUEAlgoT<T, STRATEGY>::calculateLocalDensity(const T& lt, const unsigned int layerId, float delta) {
0284 if constexpr (std::is_same_v<STRATEGY, HGCalSiliconStrategy>) {
0285 calculateLocalDensity(lt, layerId, delta, HGCalSiliconStrategy());
0286 } else {
0287 calculateLocalDensity(lt, layerId, delta, HGCalScintillatorStrategy());
0288 }
0289 }
0290
0291 template <typename T, typename STRATEGY>
0292 void HGCalCLUEAlgoT<T, STRATEGY>::calculateDistanceToHigher(const T& lt, const unsigned int layerId, float delta) {
0293 auto& cellsOnLayer = cells_[layerId];
0294 unsigned int numberOfCells = cellsOnLayer.detid.size();
0295
0296 for (unsigned int i = 0; i < numberOfCells; i++) {
0297
0298 float maxDelta = std::numeric_limits<float>::max();
0299 float i_delta = maxDelta;
0300 int i_nearestHigher = -1;
0301 auto range = outlierDeltaFactor_ * delta;
0302 std::array<int, 4> search_box = lt.searchBox(cellsOnLayer.dim1[i] - range,
0303 cellsOnLayer.dim1[i] + range,
0304 cellsOnLayer.dim2[i] - range,
0305 cellsOnLayer.dim2[i] + range);
0306
0307 for (int dim1Bin = search_box[0]; dim1Bin < search_box[1] + 1; ++dim1Bin) {
0308 for (int dim2Bin = search_box[2]; dim2Bin < search_box[3] + 1; ++dim2Bin) {
0309
0310 size_t binId = lt.getGlobalBinByBin(dim1Bin, dim2Bin);
0311 if constexpr (std::is_same_v<STRATEGY, HGCalScintillatorStrategy>)
0312 binId = lt.getGlobalBinByBin(dim1Bin, (dim2Bin % T::type::nRows));
0313
0314 size_t binSize = lt[binId].size();
0315
0316
0317 for (unsigned int j = 0; j < binSize; j++) {
0318 unsigned int otherId = lt[binId][j];
0319 float dist = distance(lt, i, otherId, layerId);
0320 bool foundHigher =
0321 (cellsOnLayer.rho[otherId] > cellsOnLayer.rho[i]) ||
0322 (cellsOnLayer.rho[otherId] == cellsOnLayer.rho[i] && cellsOnLayer.detid[otherId] > cellsOnLayer.detid[i]);
0323 if (foundHigher && dist <= i_delta) {
0324
0325 i_delta = dist;
0326
0327 i_nearestHigher = otherId;
0328 }
0329 }
0330 }
0331 }
0332 bool foundNearestHigherInSearchBox = (i_delta != maxDelta);
0333 if (foundNearestHigherInSearchBox) {
0334 cellsOnLayer.delta[i] = i_delta;
0335 cellsOnLayer.nearestHigher[i] = i_nearestHigher;
0336 } else {
0337
0338
0339 cellsOnLayer.delta[i] = maxDelta;
0340 cellsOnLayer.nearestHigher[i] = -1;
0341 }
0342
0343 LogDebug("HGCalCLUEAlgo") << "Debugging calculateDistanceToHigher: \n"
0344 << " cell: " << i << " eta: " << cellsOnLayer.dim1[i] << " phi: " << cellsOnLayer.dim2[i]
0345 << " energy: " << cellsOnLayer.weight[i] << " density: " << cellsOnLayer.rho[i]
0346 << " nearest higher: " << cellsOnLayer.nearestHigher[i]
0347 << " distance: " << cellsOnLayer.delta[i] << "\n";
0348 }
0349 }
0350
0351 template <typename T, typename STRATEGY>
0352 int HGCalCLUEAlgoT<T, STRATEGY>::findAndAssignClusters(const unsigned int layerId, float delta) {
0353
0354
0355
0356
0357 unsigned int nClustersOnLayer = 0;
0358 auto& cellsOnLayer = cells_[layerId];
0359 unsigned int numberOfCells = cellsOnLayer.detid.size();
0360 std::vector<int> localStack;
0361
0362 for (unsigned int i = 0; i < numberOfCells; i++) {
0363 float rho_c = kappa_ * cellsOnLayer.sigmaNoise[i];
0364
0365 cellsOnLayer.clusterIndex[i] = -1;
0366 bool isSeed = (cellsOnLayer.delta[i] > delta) && (cellsOnLayer.rho[i] >= rho_c);
0367 bool isOutlier = (cellsOnLayer.delta[i] > outlierDeltaFactor_ * delta) && (cellsOnLayer.rho[i] < rho_c);
0368 if (isSeed) {
0369 cellsOnLayer.clusterIndex[i] = nClustersOnLayer;
0370 cellsOnLayer.isSeed[i] = true;
0371 nClustersOnLayer++;
0372 localStack.push_back(i);
0373
0374 } else if (!isOutlier) {
0375 cellsOnLayer.followers[cellsOnLayer.nearestHigher[i]].push_back(i);
0376 }
0377 }
0378
0379
0380 while (!localStack.empty()) {
0381 int endStack = localStack.back();
0382 auto& thisSeed = cellsOnLayer.followers[endStack];
0383 localStack.pop_back();
0384
0385
0386 for (int j : thisSeed) {
0387
0388 cellsOnLayer.clusterIndex[j] = cellsOnLayer.clusterIndex[endStack];
0389
0390 localStack.push_back(j);
0391 }
0392 }
0393 return nClustersOnLayer;
0394 }
0395
0396 template <typename T, typename STRATEGY>
0397 void HGCalCLUEAlgoT<T, STRATEGY>::computeThreshold() {
0398
0399
0400
0401
0402
0403
0404
0405
0406 if (initialized_)
0407 return;
0408
0409 initialized_ = true;
0410
0411 std::vector<double> dummy;
0412
0413 dummy.resize(maxNumberOfThickIndices_ + !isNose_, 0);
0414 thresholds_.resize(maxlayer_, dummy);
0415 v_sigmaNoise_.resize(maxlayer_, dummy);
0416
0417 for (unsigned ilayer = 1; ilayer <= maxlayer_; ++ilayer) {
0418 for (unsigned ithick = 0; ithick < maxNumberOfThickIndices_; ++ithick) {
0419 float sigmaNoise = 0.001f * fcPerEle_ * nonAgedNoises_[ithick] * dEdXweights_[ilayer] /
0420 (fcPerMip_[ithick] * thicknessCorrection_[ithick]);
0421 thresholds_[ilayer - 1][ithick] = sigmaNoise * ecut_;
0422 v_sigmaNoise_[ilayer - 1][ithick] = sigmaNoise;
0423 LogDebug("HGCalCLUEAlgo") << "ilayer: " << ilayer << " nonAgedNoises: " << nonAgedNoises_[ithick]
0424 << " fcPerEle: " << fcPerEle_ << " fcPerMip: " << fcPerMip_[ithick]
0425 << " noiseMip: " << fcPerEle_ * nonAgedNoises_[ithick] / fcPerMip_[ithick]
0426 << " sigmaNoise: " << sigmaNoise << "\n";
0427 }
0428
0429 if (!isNose_) {
0430 float scintillators_sigmaNoise = 0.001f * noiseMip_ * dEdXweights_[ilayer] / sciThicknessCorrection_;
0431 thresholds_[ilayer - 1][maxNumberOfThickIndices_] = ecut_ * scintillators_sigmaNoise;
0432 v_sigmaNoise_[ilayer - 1][maxNumberOfThickIndices_] = scintillators_sigmaNoise;
0433 LogDebug("HGCalCLUEAlgo") << "ilayer: " << ilayer << " noiseMip: " << noiseMip_
0434 << " scintillators_sigmaNoise: " << scintillators_sigmaNoise << "\n";
0435 }
0436 }
0437 }
0438
0439
0440 template class HGCalCLUEAlgoT<HGCalSiliconLayerTiles, HGCalSiliconStrategy>;
0441 template class HGCalCLUEAlgoT<HGCalScintillatorLayerTiles, HGCalScintillatorStrategy>;
0442 template class HGCalCLUEAlgoT<HFNoseLayerTiles, HGCalSiliconStrategy>;