src/backend/HGCalHistoClusteringImpl.cc

0001 #include "L1Trigger/L1THGCal/interface/backend/HGCalHistoClusteringImpl.h"
0002 #include "L1Trigger/L1THGCal/interface/backend/HGCalShowerShape.h"
0003 #include "DataFormats/Math/interface/deltaR.h"
0004 #include "DataFormats/Math/interface/deltaPhi.h"
0005
0006 HGCalHistoClusteringImpl::HGCalHistoClusteringImpl(const edm::ParameterSet& conf)
0007     : dr_(conf.getParameter<double>("dR_multicluster")),
0008       dr_byLayer_coefficientA_(conf.existsAs<std::vector<double>>("dR_multicluster_byLayer_coefficientA")
0009                                    ? conf.getParameter<std::vector<double>>("dR_multicluster_byLayer_coefficientA")
0010                                    : std::vector<double>()),
0011       dr_byLayer_coefficientB_(conf.existsAs<std::vector<double>>("dR_multicluster_byLayer_coefficientB")
0012                                    ? conf.getParameter<std::vector<double>>("dR_multicluster_byLayer_coefficientB")
0013                                    : std::vector<double>()),
0014       ptC3dThreshold_(conf.getParameter<double>("minPt_multicluster")),
0015       cluster_association_input_(conf.getParameter<string>("cluster_association")),
0016       shape_(conf) {
0017   if (cluster_association_input_ == "NearestNeighbour") {
0018     cluster_association_strategy_ = NearestNeighbour;
0019   } else if (cluster_association_input_ == "EnergySplit") {
0020     cluster_association_strategy_ = EnergySplit;
0021   } else {
0022     throw cms::Exception("HGCTriggerParameterError")
0023         << "Unknown cluster association strategy'" << cluster_association_strategy_;
0024   }
0025
0026   edm::LogInfo("HGCalMulticlusterParameters")
0027       << "Multicluster dR: " << dr_ << "\nMulticluster minimum transverse-momentum: " << ptC3dThreshold_;
0028
0029   id_ = std::unique_ptr<HGCalTriggerClusterIdentificationBase>{
0030       HGCalTriggerClusterIdentificationFactory::get()->create("HGCalTriggerClusterIdentificationBDT")};
0031   id_->initialize(conf.getParameter<edm::ParameterSet>("EGIdentification"));
0032 }
0033
0034 float HGCalHistoClusteringImpl::dR(const l1t::HGCalCluster& clu, const GlobalPoint& seed) const {
0035   return (seed - clu.centreProj()).mag();
0036 }
0037
0038 std::vector<l1t::HGCalMulticluster> HGCalHistoClusteringImpl::clusterSeedMulticluster(
0039     const std::vector<edm::Ptr<l1t::HGCalCluster>>& clustersPtrs,
0040     const std::vector<std::pair<GlobalPoint, double>>& seeds,
0041     std::vector<l1t::HGCalCluster>& rejected_clusters) const {
0042   std::map<int, l1t::HGCalMulticluster> mapSeedMulticluster;
0043   std::vector<l1t::HGCalMulticluster> multiclustersOut;
0044
0045   for (const auto& clu : clustersPtrs) {
0046     int z_side = triggerTools_.zside(clu->detId());
0047
0048     double radiusCoefficientA =
0049         dr_byLayer_coefficientA_.empty() ? dr_ : dr_byLayer_coefficientA_[triggerTools_.layerWithOffset(clu->detId())];
0050     double radiusCoefficientB =
0051         dr_byLayer_coefficientB_.empty() ? 0 : dr_byLayer_coefficientB_[triggerTools_.layerWithOffset(clu->detId())];
0052
0053     double minDist = radiusCoefficientA + radiusCoefficientB * (kMidRadius_ - std::abs(clu->eta()));
0054
0055     std::vector<pair<int, double>> targetSeedsEnergy;
0056
0057     for (unsigned int iseed = 0; iseed < seeds.size(); iseed++) {
0058       GlobalPoint seedPosition = seeds[iseed].first;
0059       double seedEnergy = seeds[iseed].second;
0060
0061       if (z_side * seedPosition.z() < 0)
0062         continue;
0063       double d = this->dR(*clu, seeds[iseed].first);
0064
0065       if (d < minDist) {
0066         if (cluster_association_strategy_ == EnergySplit) {
0067           targetSeedsEnergy.emplace_back(iseed, seedEnergy);
0068         } else if (cluster_association_strategy_ == NearestNeighbour) {
0069           minDist = d;
0070
0071           if (targetSeedsEnergy.empty()) {
0072             targetSeedsEnergy.emplace_back(iseed, seedEnergy);
0073           } else {
0074             targetSeedsEnergy.at(0).first = iseed;
0075             targetSeedsEnergy.at(0).second = seedEnergy;
0076           }
0077         }
0078       }
0079     }
0080
0081     if (targetSeedsEnergy.empty()) {
0082       rejected_clusters.emplace_back(*clu);
0083       continue;
0084     }
0085     //Loop over target seeds and divide up the clusters energy
0086     double totalTargetSeedEnergy = 0;
0087     for (const auto& energy : targetSeedsEnergy) {
0088       totalTargetSeedEnergy += energy.second;
0089     }
0090
0091     for (const auto& energy : targetSeedsEnergy) {
0092       double seedWeight = 1;
0093       if (cluster_association_strategy_ == EnergySplit)
0094         seedWeight = energy.second / totalTargetSeedEnergy;
0095       if (mapSeedMulticluster[energy.first].size() == 0) {
0096         mapSeedMulticluster[energy.first] = l1t::HGCalMulticluster(clu, seedWeight);
0097       } else {
0098         mapSeedMulticluster[energy.first].addConstituent(clu, true, seedWeight);
0099       }
0100     }
0101   }
0102
0103   multiclustersOut.reserve(mapSeedMulticluster.size());
0104   for (const auto& mclu : mapSeedMulticluster)
0105     multiclustersOut.emplace_back(mclu.second);
0106
0107   return multiclustersOut;
0108 }
0109
0110 void HGCalHistoClusteringImpl::clusterizeHisto(const std::vector<edm::Ptr<l1t::HGCalCluster>>& clustersPtrs,
0111                                                const std::vector<std::pair<GlobalPoint, double>>& seedPositionsEnergy,
0112                                                const HGCalTriggerGeometryBase& triggerGeometry,
0113                                                l1t::HGCalMulticlusterBxCollection& multiclusters,
0114                                                l1t::HGCalClusterBxCollection& rejected_clusters) const {
0115   /* clusterize clusters around seeds */
0116   std::vector<l1t::HGCalCluster> rejected_clusters_vec;
0117   std::vector<l1t::HGCalMulticluster> multiclusters_vec =
0118       clusterSeedMulticluster(clustersPtrs, seedPositionsEnergy, rejected_clusters_vec);
0119   /* making the collection of multiclusters */
0120   finalizeClusters(multiclusters_vec, rejected_clusters_vec, multiclusters, rejected_clusters, triggerGeometry);
0121 }
0122
0123 void HGCalHistoClusteringImpl::finalizeClusters(std::vector<l1t::HGCalMulticluster>& multiclusters_in,
0124                                                 const std::vector<l1t::HGCalCluster>& rejected_clusters_in,
0125                                                 l1t::HGCalMulticlusterBxCollection& multiclusters_out,
0126                                                 l1t::HGCalClusterBxCollection& rejected_clusters_out,
0127                                                 const HGCalTriggerGeometryBase& triggerGeometry) const {
0128   for (const auto& tc : rejected_clusters_in) {
0129     rejected_clusters_out.push_back(0, tc);
0130   }
0131
0132   for (auto& multicluster : multiclusters_in) {
0133     // compute the eta, phi from its barycenter
0134     // + pT as scalar sum of pT of constituents
0135     double sumPt = multicluster.sumPt();
0136
0137     math::PtEtaPhiMLorentzVector multiclusterP4(sumPt, multicluster.centre().eta(), multicluster.centre().phi(), 0.);
0138     multicluster.setP4(multiclusterP4);
0139
0140     if (multicluster.pt() > ptC3dThreshold_) {
0141       //compute shower shapes
0142       shape_.fillShapes(multicluster, triggerGeometry);
0143       // fill quality flag
0144       unsigned hwQual = 0;
0145       for (unsigned wp = 0; wp < id_->working_points().size(); wp++) {
0146         hwQual |= (id_->decision(multicluster, wp) << wp);
0147       }
0148       multicluster.setHwQual(hwQual);
0149       // fill H/E
0150       multicluster.saveHOverE();
0151
0152       multiclusters_out.push_back(0, multicluster);
0153     } else {
0154       for (const auto& tc : multicluster.constituents()) {
0155         rejected_clusters_out.push_back(0, *(tc.second));
0156       }
0157     }
0158   }
0159 }