File indexing completed on 2023-03-17 11:12:21
0001 #include <unordered_set>
0002 #include <unordered_map>
0003 #include "L1Trigger/L1THGCal/interface/backend/HGCalClusteringImpl.h"
0004 #include "DataFormats/Common/interface/PtrVector.h"
0005 #include "DataFormats/Common/interface/OrphanHandle.h"
0006
0007
0008 HGCalClusteringImpl::HGCalClusteringImpl(const edm::ParameterSet& conf)
0009 : siliconSeedThreshold_(conf.getParameter<double>("seeding_threshold_silicon")),
0010 siliconTriggerCellThreshold_(conf.getParameter<double>("clustering_threshold_silicon")),
0011 scintillatorSeedThreshold_(conf.getParameter<double>("seeding_threshold_scintillator")),
0012 scintillatorTriggerCellThreshold_(conf.getParameter<double>("clustering_threshold_scintillator")),
0013 dr_(conf.getParameter<double>("dR_cluster")),
0014 clusteringAlgorithmType_(conf.getParameter<string>("clusterType")),
0015 calibSF_(conf.getParameter<double>("calibSF_cluster")),
0016 layerWeights_(conf.getParameter<std::vector<double>>("layerWeights")),
0017 applyLayerWeights_(conf.getParameter<bool>("applyLayerCalibration")) {
0018 edm::LogInfo("HGCalClusterParameters") << "C2d Clustering Algorithm selected : " << clusteringAlgorithmType_;
0019 edm::LogInfo("HGCalClusterParameters") << "C2d silicon seeding Thr: " << siliconSeedThreshold_;
0020 edm::LogInfo("HGCalClusterParameters") << "C2d silicon clustering Thr: " << siliconTriggerCellThreshold_;
0021 edm::LogInfo("HGCalClusterParameters") << "C2d scintillator seeding Thr: " << scintillatorSeedThreshold_;
0022 edm::LogInfo("HGCalClusterParameters") << "C2d scintillator clustering Thr: " << scintillatorTriggerCellThreshold_;
0023 edm::LogInfo("HGCalClusterParameters") << "C2d global calibration factor: " << calibSF_;
0024 }
0025
0026
0027 bool HGCalClusteringImpl::isPertinent(const l1t::HGCalTriggerCell& tc,
0028 const l1t::HGCalCluster& clu,
0029 double distXY) const {
0030 DetId tcDetId(tc.detId());
0031 DetId cluDetId(clu.detId());
0032 if ((triggerTools_.layer(tcDetId) != triggerTools_.layer(cluDetId)) || (tcDetId.subdetId() != cluDetId.subdetId()) ||
0033 (triggerTools_.zside(tcDetId) != triggerTools_.zside(cluDetId))) {
0034 return false;
0035 }
0036 if (clu.distance((tc)) < distXY) {
0037 return true;
0038 }
0039 return false;
0040 }
0041
0042 void HGCalClusteringImpl::clusterizeDR(const std::vector<edm::Ptr<l1t::HGCalTriggerCell>>& triggerCellsPtrs,
0043 l1t::HGCalClusterBxCollection& clusters) {
0044 bool isSeed[triggerCellsPtrs.size()];
0045
0046
0047 int itc(0);
0048 for (std::vector<edm::Ptr<l1t::HGCalTriggerCell>>::const_iterator tc = triggerCellsPtrs.begin();
0049 tc != triggerCellsPtrs.end();
0050 ++tc, ++itc) {
0051 double seedThreshold = ((*tc)->subdetId() == HGCHEB ? scintillatorSeedThreshold_ : siliconSeedThreshold_);
0052 isSeed[itc] = ((*tc)->mipPt() > seedThreshold) ? true : false;
0053 }
0054
0055
0056 std::vector<l1t::HGCalCluster> clustersTmp;
0057
0058 itc = 0;
0059 for (std::vector<edm::Ptr<l1t::HGCalTriggerCell>>::const_iterator tc = triggerCellsPtrs.begin();
0060 tc != triggerCellsPtrs.end();
0061 ++tc, ++itc) {
0062 double threshold = ((*tc)->subdetId() == HGCHEB ? scintillatorTriggerCellThreshold_ : siliconTriggerCellThreshold_);
0063 if ((*tc)->mipPt() < threshold) {
0064 continue;
0065 }
0066
0067
0068
0069
0070
0071
0072
0073 double minDist = dr_;
0074 int targetClu = -1;
0075
0076 for (unsigned iclu = 0; iclu < clustersTmp.size(); iclu++) {
0077 if (!this->isPertinent(**tc, clustersTmp.at(iclu), dr_))
0078 continue;
0079
0080 double d = clustersTmp.at(iclu).distance(**tc);
0081 if (d < minDist) {
0082 minDist = d;
0083 targetClu = int(iclu);
0084 }
0085 }
0086
0087 if (targetClu < 0 && isSeed[itc])
0088 clustersTmp.emplace_back(*tc);
0089 else if (targetClu >= 0)
0090 clustersTmp.at(targetClu).addConstituent(*tc);
0091 }
0092
0093
0094 clusters.resize(0, clustersTmp.size());
0095 for (unsigned i(0); i < clustersTmp.size(); ++i) {
0096 calibratePt(clustersTmp.at(i));
0097 clusters.set(0, i, clustersTmp.at(i));
0098 }
0099 }
0100
0101
0102
0103
0104 void HGCalClusteringImpl::triggerCellReshuffling(
0105 const std::vector<edm::Ptr<l1t::HGCalTriggerCell>>& triggerCellsPtrs,
0106 std::array<std::vector<std::vector<edm::Ptr<l1t::HGCalTriggerCell>>>, kNSides_>& reshuffledTriggerCells) {
0107 for (const auto& tc : triggerCellsPtrs) {
0108 DetId tcDetId(tc->detId());
0109 int endcap = triggerTools_.zside(tcDetId) == -1 ? 0 : 1;
0110 unsigned layer = triggerTools_.layerWithOffset(tc->detId());
0111
0112 reshuffledTriggerCells[endcap][layer - 1].emplace_back(tc);
0113 }
0114 }
0115
0116
0117 void HGCalClusteringImpl::mergeClusters(l1t::HGCalCluster& main_cluster,
0118 const l1t::HGCalCluster& secondary_cluster) const {
0119 const std::unordered_map<uint32_t, edm::Ptr<l1t::HGCalTriggerCell>>& pertinentTC = secondary_cluster.constituents();
0120
0121 for (const auto& id_tc : pertinentTC) {
0122 main_cluster.addConstituent(id_tc.second);
0123 }
0124 }
0125
0126 void HGCalClusteringImpl::NNKernel(const std::vector<edm::Ptr<l1t::HGCalTriggerCell>>& reshuffledTriggerCells,
0127 l1t::HGCalClusterBxCollection& clusters,
0128 const HGCalTriggerGeometryBase& triggerGeometry) {
0129
0130 std::vector<l1t::HGCalCluster> clustersTmp;
0131
0132
0133 std::unordered_map<uint32_t, unsigned> cluNNmap;
0134
0135
0136 for (const auto& tc_ptr : reshuffledTriggerCells) {
0137 double threshold =
0138 (tc_ptr->subdetId() == HGCHEB ? scintillatorTriggerCellThreshold_ : siliconTriggerCellThreshold_);
0139 if (tc_ptr->mipPt() < threshold) {
0140 continue;
0141 }
0142
0143
0144
0145
0146 bool createNewC2d(true);
0147 const auto neighbors = triggerGeometry.getNeighborsFromTriggerCell(tc_ptr->detId());
0148 for (const auto neighbor : neighbors) {
0149 auto tc_cluster_itr = cluNNmap.find(neighbor);
0150 if (tc_cluster_itr != cluNNmap.end()) {
0151 createNewC2d = false;
0152 if (tc_cluster_itr->second < clustersTmp.size()) {
0153 clustersTmp.at(tc_cluster_itr->second).addConstituent(tc_ptr);
0154
0155 cluNNmap.emplace(tc_ptr->detId(), tc_cluster_itr->second);
0156 } else {
0157 throw cms::Exception("HGCTriggerUnexpected")
0158 << "Trying to access a non-existing cluster. But it should exist...\n";
0159 }
0160 break;
0161 }
0162 }
0163 if (createNewC2d) {
0164 clustersTmp.emplace_back(tc_ptr);
0165 clustersTmp.back().setValid(true);
0166
0167 cluNNmap.emplace(tc_ptr->detId(), clustersTmp.size() - 1);
0168 }
0169 }
0170
0171
0172
0173 for (auto& cluster1 : clustersTmp) {
0174
0175 if (!cluster1.valid())
0176 continue;
0177
0178
0179 std::unordered_set<uint32_t> cluTcSet;
0180 for (const auto& tc_clu1 : cluster1.constituents()) {
0181 cluTcSet.insert(tc_clu1.second->detId());
0182 const auto neighbors = triggerGeometry.getNeighborsFromTriggerCell(tc_clu1.second->detId());
0183 for (const auto neighbor : neighbors) {
0184 cluTcSet.insert(neighbor);
0185 }
0186 }
0187
0188 for (auto& cluster2 : clustersTmp) {
0189
0190 if (cluster1.detId() == cluster2.detId())
0191 continue;
0192 if (!cluster2.valid())
0193 continue;
0194
0195
0196 for (const auto& tc_clu2 : cluster2.constituents()) {
0197 if (cluTcSet.find(tc_clu2.second->detId()) != cluTcSet.end()) {
0198 mergeClusters(cluster1, cluster2);
0199 cluTcSet.insert(tc_clu2.second->detId());
0200 const auto neighbors = triggerGeometry.getNeighborsFromTriggerCell(tc_clu2.second->detId());
0201 for (const auto neighbor : neighbors) {
0202 cluTcSet.insert(neighbor);
0203 }
0204 cluster2.setValid(false);
0205 break;
0206 }
0207 }
0208 }
0209 }
0210
0211
0212
0213 for (auto& cluster : clustersTmp) {
0214 if (!cluster.valid())
0215 continue;
0216 bool saveInCollection(false);
0217 for (const auto& id_tc : cluster.constituents()) {
0218
0219 double seedThreshold = (id_tc.second->subdetId() == HGCHEB ? scintillatorSeedThreshold_ : siliconSeedThreshold_);
0220 if (id_tc.second->mipPt() > seedThreshold) {
0221 saveInCollection = true;
0222 break;
0223 }
0224 }
0225 if (saveInCollection) {
0226 calibratePt(cluster);
0227 clusters.push_back(0, cluster);
0228 }
0229 }
0230 }
0231
0232 void HGCalClusteringImpl::clusterizeNN(const std::vector<edm::Ptr<l1t::HGCalTriggerCell>>& triggerCellsPtrs,
0233 l1t::HGCalClusterBxCollection& clusters,
0234 const HGCalTriggerGeometryBase& triggerGeometry) {
0235 std::array<std::vector<std::vector<edm::Ptr<l1t::HGCalTriggerCell>>>, kNSides_> reshuffledTriggerCells;
0236 unsigned layers = triggerTools_.layers(ForwardSubdetector::ForwardEmpty);
0237 for (unsigned side = 0; side < kNSides_; side++) {
0238 reshuffledTriggerCells[side].resize(layers);
0239 }
0240 triggerCellReshuffling(triggerCellsPtrs, reshuffledTriggerCells);
0241
0242 for (unsigned iec = 0; iec < kNSides_; ++iec) {
0243 for (unsigned il = 0; il < layers; ++il) {
0244 NNKernel(reshuffledTriggerCells[iec][il], clusters, triggerGeometry);
0245 }
0246 }
0247 }
0248
0249
0250 void HGCalClusteringImpl::clusterizeDRNN(const std::vector<edm::Ptr<l1t::HGCalTriggerCell>>& triggerCellsPtrs,
0251 l1t::HGCalClusterBxCollection& clusters,
0252 const HGCalTriggerGeometryBase& triggerGeometry) {
0253 bool isSeed[triggerCellsPtrs.size()];
0254 std::vector<unsigned> seedPositions;
0255 seedPositions.reserve(triggerCellsPtrs.size());
0256
0257
0258 int itc(0);
0259 for (std::vector<edm::Ptr<l1t::HGCalTriggerCell>>::const_iterator tc = triggerCellsPtrs.begin();
0260 tc != triggerCellsPtrs.end();
0261 ++tc, ++itc) {
0262 double seedThreshold = ((*tc)->subdetId() == HGCHEB ? scintillatorSeedThreshold_ : siliconSeedThreshold_);
0263
0264
0265 isSeed[itc] = ((*tc)->mipPt() > seedThreshold) ? true : false;
0266 if (isSeed[itc]) {
0267 seedPositions.push_back(itc);
0268
0269
0270 for (auto pos : seedPositions) {
0271 if (((*tc)->position() - triggerCellsPtrs[pos]->position()).mag() < dr_) {
0272 if (this->areTCneighbour((*tc)->detId(), triggerCellsPtrs[pos]->detId(), triggerGeometry)) {
0273 isSeed[itc] = false;
0274 seedPositions.pop_back();
0275 }
0276 }
0277 }
0278 }
0279 }
0280
0281
0282 std::vector<l1t::HGCalCluster> clustersTmp;
0283
0284
0285 clustersTmp.reserve(seedPositions.size());
0286 for (auto pos : seedPositions) {
0287 clustersTmp.emplace_back(triggerCellsPtrs[pos]);
0288 }
0289
0290
0291 itc = 0;
0292 for (std::vector<edm::Ptr<l1t::HGCalTriggerCell>>::const_iterator tc = triggerCellsPtrs.begin();
0293 tc != triggerCellsPtrs.end();
0294 ++tc, ++itc) {
0295
0296 double threshold = ((*tc)->subdetId() == HGCHEB ? scintillatorTriggerCellThreshold_ : siliconTriggerCellThreshold_);
0297
0298
0299 if ((*tc)->mipPt() < threshold)
0300 continue;
0301 if (isSeed[itc])
0302 continue;
0303
0304
0305 std::vector<unsigned> tcPertinentClusters;
0306 unsigned iclu(0);
0307
0308 for (const auto& clu : clustersTmp) {
0309 if (this->isPertinent(**tc, clu, dr_)) {
0310 tcPertinentClusters.push_back(iclu);
0311 }
0312 ++iclu;
0313 }
0314
0315 if (tcPertinentClusters.empty()) {
0316 continue;
0317 } else if (tcPertinentClusters.size() == 1) {
0318 clustersTmp.at(tcPertinentClusters.at(0)).addConstituent(*tc);
0319 } else {
0320
0321 double totMipt = 0;
0322 for (auto clu : tcPertinentClusters) {
0323 totMipt += clustersTmp.at(clu).seedMipPt();
0324 }
0325
0326 for (auto clu : tcPertinentClusters) {
0327 double seedMipt = clustersTmp.at(clu).seedMipPt();
0328 clustersTmp.at(clu).addConstituent(*tc, true, seedMipt / totMipt);
0329 }
0330 }
0331 }
0332
0333
0334 clusters.resize(0, clustersTmp.size());
0335 for (unsigned i(0); i < clustersTmp.size(); ++i) {
0336 this->removeUnconnectedTCinCluster(clustersTmp.at(i), triggerGeometry);
0337 calibratePt(clustersTmp.at(i));
0338 clusters.set(0, i, clustersTmp.at(i));
0339 }
0340 }
0341
0342 bool HGCalClusteringImpl::areTCneighbour(uint32_t detIDa,
0343 uint32_t detIDb,
0344 const HGCalTriggerGeometryBase& triggerGeometry) {
0345 const auto neighbors = triggerGeometry.getNeighborsFromTriggerCell(detIDa);
0346
0347 if (neighbors.find(detIDb) != neighbors.end())
0348 return true;
0349
0350 return false;
0351 }
0352
0353 void HGCalClusteringImpl::removeUnconnectedTCinCluster(l1t::HGCalCluster& cluster,
0354 const HGCalTriggerGeometryBase& triggerGeometry) {
0355
0356 const std::unordered_map<uint32_t, edm::Ptr<l1t::HGCalTriggerCell>>& constituents = cluster.constituents();
0357 Basic3DVector<float> seedCentre(constituents.at(cluster.detId())->position());
0358
0359
0360 vector<pair<edm::Ptr<l1t::HGCalTriggerCell>, float>> distances;
0361 for (const auto& id_tc : constituents) {
0362 Basic3DVector<float> tcCentre(id_tc.second->position());
0363 float distance = (seedCentre - tcCentre).mag();
0364 distances.push_back(pair<edm::Ptr<l1t::HGCalTriggerCell>, float>(id_tc.second, distance));
0365 }
0366
0367
0368
0369 std::sort(distances.begin(), distances.end(), distanceSorter);
0370
0371
0372 bool toRemove[constituents.size()];
0373 toRemove[0] = false;
0374 for (unsigned itc = 1; itc < distances.size(); itc++) {
0375
0376 toRemove[itc] = true;
0377 const edm::Ptr<l1t::HGCalTriggerCell>& tcToStudy = distances[itc].first;
0378
0379
0380 for (unsigned itc_ref = 0; itc_ref < itc; itc_ref++) {
0381 if (!toRemove[itc_ref]) {
0382 if (areTCneighbour(tcToStudy->detId(), distances.at(itc_ref).first->detId(), triggerGeometry)) {
0383 toRemove[itc] = false;
0384 break;
0385 }
0386 }
0387 }
0388 }
0389
0390
0391 for (unsigned i = 0; i < distances.size(); i++) {
0392 if (toRemove[i])
0393 cluster.removeConstituent(distances.at(i).first);
0394 }
0395 }
0396
0397 void HGCalClusteringImpl::calibratePt(l1t::HGCalCluster& cluster) {
0398 double calibPt = 0.;
0399
0400 if (applyLayerWeights_) {
0401 unsigned layerN = triggerTools_.layerWithOffset(cluster.detId());
0402
0403 if (layerWeights_.at(layerN) == 0.) {
0404 throw cms::Exception("BadConfiguration")
0405 << "2D cluster energy forced to 0 by calibration coefficients.\n"
0406 << "The configuration should be changed. "
0407 << "Discarded layers should be defined in hgcalTriggerGeometryESProducer.TriggerGeometry.DisconnectedLayers "
0408 "and not with calibration coefficients = 0\n";
0409 }
0410
0411 calibPt = layerWeights_.at(layerN) * cluster.mipPt();
0412
0413 }
0414
0415 else {
0416 calibPt = cluster.pt() * calibSF_;
0417 }
0418
0419 math::PtEtaPhiMLorentzVector calibP4(calibPt, cluster.eta(), cluster.phi(), 0.);
0420
0421 cluster.setP4(calibP4);
0422 }