Back to home page

Project CMSSW displayed by LXR
 
 

    

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 //class constructor
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 /* dR-algorithms */
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   /* search for cluster seeds */
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   /* clustering the TCs */
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            searching for TC near the center of the cluster  
0069            ToBeFixed: if a tc is not a seed, but could be potencially be part of a cluster generated by a late seed, 
0070                       the tc will not be clusterized  
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   /* store clusters in the persistent collection */
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 /* NN-algorithms */
0102 
0103 /* storing trigger cells into vector per layer and per endcap */
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 /* merge clusters that have common neighbors */
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   /* declaring the clusters vector */
0130   std::vector<l1t::HGCalCluster> clustersTmp;
0131 
0132   // map TC id -> cluster index in clustersTmp
0133   std::unordered_map<uint32_t, unsigned> cluNNmap;
0134 
0135   /* loop over the trigger-cells */
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     // Check if the neighbors of that TC are already included in a cluster
0144     // If this is the case, add the TC to the first (arbitrary) neighbor cluster
0145     // Otherwise create a new cluster
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           // map TC id to the existing cluster
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       // map TC id to the cluster index (size - 1)
0167       cluNNmap.emplace(tc_ptr->detId(), clustersTmp.size() - 1);
0168     }
0169   }
0170 
0171   /* declaring the vector with possible clusters merged */
0172   // Merge neighbor clusters together
0173   for (auto& cluster1 : clustersTmp) {
0174     // If the cluster has been merged into another one, skip it
0175     if (!cluster1.valid())
0176       continue;
0177     // Fill a set containing all TC included in the clusters
0178     // as well as all neighbor TC
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       // If the cluster has been merged into another one, skip it
0190       if (cluster1.detId() == cluster2.detId())
0191         continue;
0192       if (!cluster2.valid())
0193         continue;
0194       // Check if the TC in clu2 are in clu1 or its neighbors
0195       // If yes, merge the second cluster into the first one
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   /* store clusters in the persistent collection */
0212   // only if the cluster contain a TC above the seed threshold
0213   for (auto& cluster : clustersTmp) {
0214     if (!cluster.valid())
0215       continue;
0216     bool saveInCollection(false);
0217     for (const auto& id_tc : cluster.constituents()) {
0218       /* threshold in transverse-mip */
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 /*** FW-algorithms ***/
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   /* search for cluster seeds */
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     /* decide if is a seed, if yes store the position into of triggerCellsPtrs */
0265     isSeed[itc] = ((*tc)->mipPt() > seedThreshold) ? true : false;
0266     if (isSeed[itc]) {
0267       seedPositions.push_back(itc);
0268 
0269       /* remove tc from the seed vector if is a NN of an other seed*/
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   /* clustering the TCs */
0282   std::vector<l1t::HGCalCluster> clustersTmp;
0283 
0284   // every seed generates a cluster
0285   clustersTmp.reserve(seedPositions.size());
0286   for (auto pos : seedPositions) {
0287     clustersTmp.emplace_back(triggerCellsPtrs[pos]);
0288   }
0289 
0290   /* add the tc to the clusters */
0291   itc = 0;
0292   for (std::vector<edm::Ptr<l1t::HGCalTriggerCell>>::const_iterator tc = triggerCellsPtrs.begin();
0293        tc != triggerCellsPtrs.end();
0294        ++tc, ++itc) {
0295     /* get the correct threshold for the different part of the detector */
0296     double threshold = ((*tc)->subdetId() == HGCHEB ? scintillatorTriggerCellThreshold_ : siliconTriggerCellThreshold_);
0297 
0298     /* continue if not passing the threshold */
0299     if ((*tc)->mipPt() < threshold)
0300       continue;
0301     if (isSeed[itc])
0302       continue;  //No sharing of seeds between clusters (TBC)
0303 
0304     /* searching for TC near the centre of the cluster  */
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       /* calculate the fractions */
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   /* store clusters in the persistent collection */
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   /* get the constituents and the centre of the seed tc (considered as the first of the constituents) */
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   /* distances from the seed */
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   /* sorting (needed in order to be sure that we are skipping any tc) */
0368   /* FIXME: better sorting needed!!! */
0369   std::sort(distances.begin(), distances.end(), distanceSorter);
0370 
0371   /* checking if the tc is connected to the seed */
0372   bool toRemove[constituents.size()];
0373   toRemove[0] = false;  // this is the seed
0374   for (unsigned itc = 1; itc < distances.size(); itc++) {
0375     /* get the tc under study */
0376     toRemove[itc] = true;
0377     const edm::Ptr<l1t::HGCalTriggerCell>& tcToStudy = distances[itc].first;
0378 
0379     /* compare with the tc in the cluster */
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   /* remove the unconnected TCs */
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 }
[ 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!