src/backend/HGCalMulticlusteringImpl.cc

0001 #include "L1Trigger/L1THGCal/interface/backend/HGCalMulticlusteringImpl.h"
0002 #include "L1Trigger/L1THGCal/interface/backend/HGCalShowerShape.h"
0003 #include "DataFormats/Math/interface/deltaR.h"
0004
0005 HGCalMulticlusteringImpl::HGCalMulticlusteringImpl(const edm::ParameterSet& conf)
0006     : dr_(conf.getParameter<double>("dR_multicluster")),
0007       ptC3dThreshold_(conf.getParameter<double>("minPt_multicluster")),
0008       multiclusterAlgoType_(conf.getParameter<string>("type_multicluster")),
0009       distDbscan_(conf.getParameter<double>("dist_dbscan_multicluster")),
0010       minNDbscan_(conf.getParameter<unsigned>("minN_dbscan_multicluster")) {
0011   edm::LogInfo("HGCalMulticlusterParameters") << "Multicluster dR for Near Neighbour search: " << dr_;
0012   edm::LogInfo("HGCalMulticlusterParameters") << "Multicluster minimum transverse-momentum: " << ptC3dThreshold_;
0013   edm::LogInfo("HGCalMulticlusterParameters") << "Multicluster DBSCAN Clustering distance: " << distDbscan_;
0014   edm::LogInfo("HGCalMulticlusterParameters") << "Multicluster clustering min number of subclusters: " << minNDbscan_;
0015   edm::LogInfo("HGCalMulticlusterParameters")
0016       << "Multicluster type of multiclustering algortihm: " << multiclusterAlgoType_;
0017   id_ = std::unique_ptr<HGCalTriggerClusterIdentificationBase>{
0018       HGCalTriggerClusterIdentificationFactory::get()->create("HGCalTriggerClusterIdentificationBDT")};
0019   id_->initialize(conf.getParameter<edm::ParameterSet>("EGIdentification"));
0020 }
0021
0022 bool HGCalMulticlusteringImpl::isPertinent(const l1t::HGCalCluster& clu,
0023                                            const l1t::HGCalMulticluster& mclu,
0024                                            double dR) const {
0025   DetId cluDetId(clu.detId());
0026   DetId firstClusterDetId(mclu.detId());
0027
0028   if (triggerTools_.zside(cluDetId) != triggerTools_.zside(firstClusterDetId)) {
0029     return false;
0030   }
0031   if ((mclu.centreProj() - clu.centreProj()).mag() < dR) {
0032     return true;
0033   }
0034   return false;
0035 }
0036
0037 void HGCalMulticlusteringImpl::findNeighbor(const std::vector<std::pair<unsigned int, double>>& rankedList,
0038                                             unsigned int searchInd,
0039                                             const std::vector<edm::Ptr<l1t::HGCalCluster>>& clustersPtrs,
0040                                             std::vector<unsigned int>& neighbors) {
0041   if (clustersPtrs.size() <= searchInd || clustersPtrs.size() < rankedList.size()) {
0042     throw cms::Exception("IndexOutOfBound: clustersPtrs in 'findNeighbor'");
0043   }
0044
0045   for (unsigned int ind = searchInd + 1;
0046        ind < rankedList.size() && fabs(rankedList.at(ind).second - rankedList.at(searchInd).second) < distDbscan_;
0047        ind++) {
0048     if (clustersPtrs.size() <= rankedList.at(ind).first) {
0049       throw cms::Exception("IndexOutOfBound: clustersPtrs in 'findNeighbor'");
0050
0051     } else if (((*(clustersPtrs[rankedList.at(ind).first])).centreProj() -
0052                 (*(clustersPtrs[rankedList.at(searchInd).first])).centreProj())
0053                    .mag() < distDbscan_) {
0054       neighbors.push_back(ind);
0055     }
0056   }
0057
0058   for (unsigned int ind = 0;
0059        ind < searchInd && fabs(rankedList.at(searchInd).second - rankedList.at(ind).second) < distDbscan_;
0060        ind++) {
0061     if (clustersPtrs.size() <= rankedList.at(ind).first) {
0062       throw cms::Exception("IndexOutOfBound: clustersPtrs in 'findNeighbor'");
0063
0064     } else if (((*(clustersPtrs[rankedList.at(ind).first])).centreProj() -
0065                 (*(clustersPtrs[rankedList.at(searchInd).first])).centreProj())
0066                    .mag() < distDbscan_) {
0067       neighbors.push_back(ind);
0068     }
0069   }
0070 }
0071
0072 void HGCalMulticlusteringImpl::clusterizeDR(const std::vector<edm::Ptr<l1t::HGCalCluster>>& clustersPtrs,
0073                                             l1t::HGCalMulticlusterBxCollection& multiclusters,
0074                                             const HGCalTriggerGeometryBase& triggerGeometry) {
0075   std::vector<l1t::HGCalMulticluster> multiclustersTmp;
0076
0077   int iclu = 0;
0078   for (std::vector<edm::Ptr<l1t::HGCalCluster>>::const_iterator clu = clustersPtrs.begin(); clu != clustersPtrs.end();
0079        ++clu, ++iclu) {
0080     double minDist = dr_;
0081     int targetMulticlu = -1;
0082
0083     for (unsigned imclu = 0; imclu < multiclustersTmp.size(); imclu++) {
0084       if (!this->isPertinent(**clu, multiclustersTmp.at(imclu), dr_))
0085         continue;
0086
0087       double d = (multiclustersTmp.at(imclu).centreProj() - (*clu)->centreProj()).mag();
0088       if (d < minDist) {
0089         minDist = d;
0090         targetMulticlu = int(imclu);
0091       }
0092     }
0093
0094     if (targetMulticlu < 0)
0095       multiclustersTmp.emplace_back(*clu);
0096     else
0097       multiclustersTmp.at(targetMulticlu).addConstituent(*clu);
0098   }
0099
0100   /* making the collection of multiclusters */
0101   finalizeClusters(multiclustersTmp, multiclusters, triggerGeometry);
0102 }
0103 void HGCalMulticlusteringImpl::clusterizeDBSCAN(const std::vector<edm::Ptr<l1t::HGCalCluster>>& clustersPtrs,
0104                                                 l1t::HGCalMulticlusterBxCollection& multiclusters,
0105                                                 const HGCalTriggerGeometryBase& triggerGeometry) {
0106   std::vector<l1t::HGCalMulticluster> multiclustersTmp;
0107   l1t::HGCalMulticluster mcluTmp;
0108   std::vector<bool> visited(clustersPtrs.size(), false);
0109   std::vector<bool> merged(clustersPtrs.size(), false);
0110   std::vector<std::pair<unsigned int, double>> rankedList;
0111   rankedList.reserve(clustersPtrs.size());
0112   std::vector<std::vector<unsigned int>> neighborList;
0113   neighborList.reserve(clustersPtrs.size());
0114
0115   int iclu = 0, imclu = 0, neighNo;
0116   double dist = 0.;
0117
0118   for (std::vector<edm::Ptr<l1t::HGCalCluster>>::const_iterator clu = clustersPtrs.begin(); clu != clustersPtrs.end();
0119        ++clu, ++iclu) {
0120     dist = (*clu)->centreProj().mag() * triggerTools_.zside((*clu)->detId());
0121     rankedList.push_back(std::make_pair(iclu, dist));
0122   }
0123   iclu = 0;
0124   std::sort(rankedList.begin(), rankedList.end(), [](auto& left, auto& right) { return left.second < right.second; });
0125
0126   for (const auto& cluRanked : rankedList) {
0127     std::vector<unsigned int> neighbors;
0128
0129     if (!visited.at(iclu)) {
0130       visited.at(iclu) = true;
0131       findNeighbor(rankedList, iclu, clustersPtrs, neighbors);
0132       neighborList.push_back(std::move(neighbors));
0133
0134       if (neighborList.at(iclu).size() >= minNDbscan_) {
0135         multiclustersTmp.emplace_back(clustersPtrs[cluRanked.first]);
0136         merged.at(iclu) = true;
0137         /* dynamic range loop: range-based loop syntax cannot be employed */
0138         for (unsigned int neighInd = 0; neighInd < neighborList.at(iclu).size(); neighInd++) {
0139           neighNo = neighborList.at(iclu).at(neighInd);
0140           /* This condition also ensures merging of clusters visited by other clusters but not merged. */
0141           if (!merged.at(neighNo)) {
0142             merged.at(neighNo) = true;
0143             multiclustersTmp.at(imclu).addConstituent(clustersPtrs[rankedList.at(neighNo).first]);
0144
0145             if (!visited.at(neighNo)) {
0146               visited.at(neighNo) = true;
0147               std::vector<unsigned int> secNeighbors;
0148               findNeighbor(rankedList, neighNo, clustersPtrs, secNeighbors);
0149
0150               if (secNeighbors.size() >= minNDbscan_) {
0151                 neighborList.at(iclu).insert(neighborList.at(iclu).end(), secNeighbors.begin(), secNeighbors.end());
0152               }
0153             }
0154           }
0155         }
0156         imclu++;
0157       }
0158     } else
0159       neighborList.push_back(std::move(neighbors));
0160     iclu++;
0161   }
0162   /* making the collection of multiclusters */
0163   finalizeClusters(multiclustersTmp, multiclusters, triggerGeometry);
0164 }
0165
0166 void HGCalMulticlusteringImpl::finalizeClusters(std::vector<l1t::HGCalMulticluster>& multiclusters_in,
0167                                                 l1t::HGCalMulticlusterBxCollection& multiclusters_out,
0168                                                 const HGCalTriggerGeometryBase& triggerGeometry) {
0169   for (auto& multicluster : multiclusters_in) {
0170     // compute the eta, phi from its barycenter
0171     // + pT as scalar sum of pT of constituents
0172     double sumPt = 0.;
0173     const std::unordered_map<uint32_t, edm::Ptr<l1t::HGCalCluster>>& clusters = multicluster.constituents();
0174     for (const auto& id_cluster : clusters)
0175       sumPt += id_cluster.second->pt();
0176
0177     math::PtEtaPhiMLorentzVector multiclusterP4(sumPt, multicluster.centre().eta(), multicluster.centre().phi(), 0.);
0178     multicluster.setP4(multiclusterP4);
0179
0180     if (multicluster.pt() > ptC3dThreshold_) {
0181       //compute shower shapes
0182       shape_.fillShapes(multicluster, triggerGeometry);
0183       // fill quality flag
0184       multicluster.setHwQual(id_->decision(multicluster));
0185       // fill H/E
0186       multicluster.saveHOverE();
0187
0188       multiclusters_out.push_back(0, multicluster);
0189     }
0190   }
0191 }