File indexing completed on 2023-03-17 11:12:21
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
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
0116 std::vector<l1t::HGCalCluster> rejected_clusters_vec;
0117 std::vector<l1t::HGCalMulticluster> multiclusters_vec =
0118 clusterSeedMulticluster(clustersPtrs, seedPositionsEnergy, rejected_clusters_vec);
0119
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
0134
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
0142 shape_.fillShapes(multicluster, triggerGeometry);
0143
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
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 }