File indexing completed on 2025-01-27 02:50:20
0001 #include "L1Trigger/Phase2L1ParticleFlow/interface/CaloClusterer.h"
0002
0003 #include <cassert>
0004
0005 #include "DataFormats/Math/interface/deltaPhi.h"
0006 #include "DataFormats/Common/interface/RefToPtr.h"
0007 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0008 #include "FWCore/Utilities/interface/Exception.h"
0009
0010 const float l1tpf_calo::Phase1Grid::phase1_towerEtas_[l1tpf_calo::Phase1Grid::phase1_nEta_] = {
0011 0, 0.087, 0.174, 0.261, 0.348, 0.435, 0.522, 0.609, 0.696, 0.783, 0.870, 0.957, 1.044, 1.131,
0012 1.218, 1.305, 1.392, 1.479, 1.566, 1.653, 1.740, 1.830, 1.930, 2.043, 2.172, 2.322, 2.5, 2.650,
0013 2.853, 3.139, 3.314, 3.489, 3.664, 3.839, 4.013, 4.191, 4.363, 4.538, 4.716, 4.889, 5.191};
0014 const float l1tpf_calo::Phase2Grid::phase2_towerEtas_[l1tpf_calo::Phase2Grid::phase2_nEta_] = {
0015 0, 0.087, 0.174, 0.261, 0.348, 0.435, 0.522, 0.609, 0.696, 0.783, 0.870, 0.957, 1.044, 1.131, 1.218, 1.305,
0016 1.392, 1.479, 1.564, 1.648, 1.732, 1.817, 1.901, 1.986, 2.071, 2.155, 2.240, 2.324, 2.409, 2.493, 2.577, 2.662,
0017 2.747, 2.831, 2.915, 3.0, 3.139, 3.314, 3.489, 3.664, 3.839, 4.013, 4.191, 4.363, 4.538, 4.716, 4.889, 5.191};
0018
0019 l1tpf_calo::Phase1GridBase::Phase1GridBase(
0020 int nEta, int nPhi, int ietaCoarse, int ietaVeryCoarse, const float *towerEtas)
0021 : Grid(2 * ((ietaCoarse - 1) * nPhi + (ietaVeryCoarse - ietaCoarse) * (nPhi / 2) +
0022 (nEta - ietaVeryCoarse + 1) * (nPhi / 4))),
0023 nEta_(nEta),
0024 nPhi_(nPhi),
0025 ietaCoarse_(ietaCoarse),
0026 ietaVeryCoarse_(ietaVeryCoarse),
0027 towerEtas_(towerEtas),
0028 cell_map_(2 * nEta * nPhi, -1) {
0029 int icell = 0;
0030 for (int ie = -nEta_; ie <= nEta_; ++ie) {
0031 int absie = std::abs(ie);
0032 for (int iph = 1; iph <= nPhi_; ++iph) {
0033 if (!valid_ieta_iphi(ie, iph))
0034 continue;
0035 ieta_[icell] = ie;
0036 iphi_[icell] = iph;
0037 float etaLo = (absie < nEta_ ? towerEtas_[absie - 1] : towerEtas_[absie - 2]);
0038 float etaHi = (absie < nEta_ ? towerEtas_[absie] : towerEtas_[absie - 1]);
0039 eta_[icell] = (ie > 0 ? 0.5 : -0.5) * (etaLo + etaHi);
0040 etaWidth_[icell] = (etaHi - etaLo);
0041 phiWidth_[icell] = 2 * M_PI / nPhi_;
0042 if (absie >= ietaVeryCoarse_)
0043 phiWidth_[icell] *= 4;
0044 else if (absie >= ietaCoarse_)
0045 phiWidth_[icell] *= 2;
0046 phi_[icell] = (iph - 1) * 2 * M_PI / nPhi_ + 0.5 * phiWidth_[icell];
0047 if (phi_[icell] > M_PI)
0048 phi_[icell] -= 2 * M_PI;
0049 std::fill(neighbours_[icell].begin(), neighbours_[icell].end(), -1);
0050 cell_map_[(ie + nEta_) + 2 * nEta_ * (iph - 1)] = icell;
0051 icell++;
0052 }
0053 }
0054 assert(unsigned(icell) == ncells_);
0055
0056 for (icell = 0; icell < int(ncells_); ++icell) {
0057 int ie = ieta_[icell], iph = iphi_[icell];
0058 int ineigh = 0;
0059 for (int deta = -1; deta <= +1; ++deta) {
0060 for (int dphi = -1; dphi <= +1; ++dphi) {
0061 if (deta == 0 && dphi == 0)
0062 continue;
0063 neighbours_[icell][ineigh++] = imove(ie, iph, deta, dphi);
0064 }
0065 }
0066 }
0067
0068
0069
0070
0071
0072
0073
0074
0075 }
0076
0077 int l1tpf_calo::Phase1GridBase::find_cell(float eta, float phi) const {
0078 int ieta =
0079 (eta != 0) ? std::distance(towerEtas_, std::lower_bound(towerEtas_, towerEtas_ + nEta_, std::abs(eta))) : 1;
0080 if (ieta == nEta_)
0081 return -1;
0082 assert(ieta > 0 && ieta < nEta_);
0083 if (ieta > nEta_)
0084 ieta = nEta_;
0085 if (eta < 0)
0086 ieta = -ieta;
0087 phi = reco::reducePhiRange(phi);
0088 if (phi < 0)
0089 phi += 2 * M_PI;
0090 int iphi = std::floor(phi * nPhi_ / (2 * M_PI));
0091 if (phi >= 2 * M_PI)
0092 iphi = nPhi_ - 1;
0093 assert(iphi < nPhi_);
0094 if (std::abs(ieta) >= ietaVeryCoarse_)
0095 iphi -= (iphi % 4);
0096 else if (std::abs(ieta) >= ietaCoarse_)
0097 iphi -= (iphi % 2);
0098 iphi += 1;
0099
0100
0101
0102
0103
0104 assert(valid_ieta_iphi(ieta, iphi));
0105 int icell = ifind_cell(ieta, iphi);
0106 assert(icell != -1);
0107
0108
0109
0110
0111
0112
0113
0114
0115 return icell;
0116 }
0117
0118 int l1tpf_calo::Phase1GridBase::imove(int ieta, int iphi, int deta, int dphi) {
0119 int ie = ieta, iph = iphi;
0120 switch (deta) {
0121 case -1:
0122 ie = (ie == -nEta_ ? 0 : (ie == +1 ? -1 : ie - 1));
0123 break;
0124 case +1:
0125 ie = (ie == +nEta_ ? 0 : (ie == -1 ? +1 : ie + 1));
0126 break;
0127 case 0:
0128 break;
0129 default:
0130 assert(false);
0131 };
0132 if (ie == 0)
0133 return -1;
0134 switch (dphi) {
0135 case -1:
0136 iph = (iph == 1 ? nPhi_ : iph - 1);
0137 break;
0138 case +1:
0139 iph = (iph == nPhi_ ? 1 : iph + 1);
0140 break;
0141 case 0:
0142 break;
0143 default:
0144 assert(false);
0145 };
0146 if (!valid_ieta_iphi(ie, iph))
0147 return -1;
0148 int icell = ifind_cell(ie, iph);
0149 assert(!(ie == ieta && iph == iphi));
0150 assert(icell != -1);
0151 assert(icell != ifind_cell(ieta, iphi));
0152 return icell;
0153 }
0154
0155 const l1tpf_calo::Grid *l1tpf_calo::getGrid(const std::string &type) {
0156 static const Phase1Grid _phase1Grid;
0157 static const Phase2Grid _phase2Grid;
0158 if (type == "phase1")
0159 return &_phase1Grid;
0160 else if (type == "phase2")
0161 return &_phase2Grid;
0162 else
0163 throw cms::Exception("Configuration") << "Unsupported grid type '" << type << "'\n";
0164 }
0165
0166 l1tpf_calo::SingleCaloClusterer::SingleCaloClusterer(const edm::ParameterSet &pset)
0167 : grid_(getGrid(pset.getParameter<std::string>("grid"))),
0168 rawet_(*grid_),
0169 unclustered_(*grid_),
0170 eta_center_(*grid_),
0171 phi_center_(*grid_),
0172 precluster_(*grid_),
0173 clusterIndex_(*grid_),
0174 cellKey_(*grid_),
0175 preciseEtaPhi_(pset.existsAs<bool>("usePreciseEtaPhi") ? pset.getParameter<bool>("usePreciseEtaPhi") : false),
0176 etaBounds_(pset.getParameter<std::vector<double>>("etaBounds")),
0177 phiBounds_(pset.getParameter<std::vector<double>>("phiBounds")),
0178 maxClustersEtaPhi_(pset.getParameter<std::vector<unsigned int>>("maxClustersEtaPhi")),
0179 clusters_(),
0180 nullCluster_(),
0181 zsEt_(pset.getParameter<double>("zsEt")),
0182 seedEt_(pset.getParameter<double>("seedEt")),
0183 minClusterEt_(pset.getParameter<double>("minClusterEt")),
0184 minEtToGrow_(pset.existsAs<double>("minEtToGrow") ? pset.getParameter<double>("minEtToGrow") : -1),
0185 energyWeightedPosition_(pset.getParameter<bool>("energyWeightedPosition")) {
0186 std::string energyShareAlgo = pset.getParameter<std::string>("energyShareAlgo");
0187 if (energyShareAlgo == "fractions")
0188 energyShareAlgo_ = EnergyShareAlgo::Fractions;
0189 else if (energyShareAlgo == "none")
0190 energyShareAlgo_ = EnergyShareAlgo::None;
0191 else if (energyShareAlgo == "greedy")
0192 energyShareAlgo_ = EnergyShareAlgo::Greedy;
0193 else if (energyShareAlgo == "crude")
0194 energyShareAlgo_ = EnergyShareAlgo::Crude;
0195 else
0196 throw cms::Exception("Configuration") << "Unsupported energyShareAlgo '" << energyShareAlgo << "'\n";
0197 if (pset.existsAs<std::vector<int>>("neighborCells")) {
0198 neighborCells_ = pset.getParameter<std::vector<int>>(
0199 "neighborCells");
0200 } else {
0201 neighborCells_ = std::vector<int>({0, 1, 2, 3, 4, 5, 6, 7});
0202
0203
0204
0205 }
0206 if ((etaBounds_.size() - 1) * (phiBounds_.size() - 1) != maxClustersEtaPhi_.size()) {
0207 throw cms::Exception("Configuration")
0208 << "Size mismatch between eta/phi bounds and max clusters: " << (etaBounds_.size() - 1) << " x "
0209 << (phiBounds_.size() - 1) << " != " << maxClustersEtaPhi_.size() << "\n";
0210 }
0211 if (!std::is_sorted(etaBounds_.begin(), etaBounds_.end())) {
0212 throw cms::Exception("Configuration") << "etaBounds is not sorted\n";
0213 }
0214 if (!std::is_sorted(phiBounds_.begin(), phiBounds_.end())) {
0215 throw cms::Exception("Configuration") << "phiBounds is not sorted\n";
0216 }
0217 }
0218
0219 l1tpf_calo::SingleCaloClusterer::~SingleCaloClusterer() {}
0220
0221 void l1tpf_calo::SingleCaloClusterer::clear() {
0222 rawet_.zero();
0223 eta_center_.zero();
0224 phi_center_.zero();
0225 clusters_.clear();
0226 clusterIndex_.fill(-1);
0227 }
0228
0229 void l1tpf_calo::SingleCaloClusterer::run() {
0230 unsigned int i, ncells = grid_->size();
0231
0232
0233 cellKey_.fill(-1);
0234 int key = 0;
0235 for (i = 0; i < ncells; ++i) {
0236 if (rawet_[i] < zsEt_) {
0237 rawet_[i] = 0;
0238 } else {
0239 cellKey_[i] = key++;
0240 }
0241 }
0242
0243 precluster_.clear();
0244
0245
0246 for (i = 0; i < ncells; ++i) {
0247 if (rawet_[i] > seedEt_) {
0248 precluster_[i].ptLocalMax = rawet_[i];
0249
0250
0251 for (const auto &ineigh : neighborCells_) {
0252 if (ineigh >= 4)
0253 continue;
0254 if (rawet_.neigh(i, ineigh) > rawet_[i])
0255 precluster_[i].ptLocalMax = 0;
0256
0257
0258
0259
0260 }
0261 for (const auto &ineigh : neighborCells_) {
0262 if (ineigh < 4)
0263 continue;
0264 if (rawet_.neigh(i, ineigh) >= rawet_[i])
0265 precluster_[i].ptLocalMax = 0;
0266
0267
0268
0269
0270 }
0271 }
0272 }
0273
0274 for (i = 0; i < ncells; ++i) {
0275 if (precluster_[i].ptLocalMax == 0) {
0276 switch (energyShareAlgo_) {
0277 case EnergyShareAlgo::Fractions: {
0278 float tot = 0;
0279 for (const auto &ineigh : neighborCells_) {
0280 tot += precluster_.neigh(i, ineigh).ptLocalMax;
0281 }
0282 precluster_[i].ptOverNeighLocalMaxSum = tot ? rawet_[i] / tot : 0;
0283 } break;
0284 case EnergyShareAlgo::None:
0285 precluster_[i].ptOverNeighLocalMaxSum = rawet_[i];
0286 break;
0287 case EnergyShareAlgo::Greedy: {
0288 float maxet = 0;
0289 for (const auto &ineigh : neighborCells_) {
0290 maxet = std::max(maxet, precluster_.neigh(i, ineigh).ptLocalMax);
0291 }
0292 precluster_[i].ptOverNeighLocalMaxSum = maxet;
0293 } break;
0294 case EnergyShareAlgo::Crude: {
0295 int number = 0;
0296 for (const auto &ineigh : neighborCells_) {
0297 number += (precluster_.neigh(i, ineigh).ptLocalMax > 0);
0298 }
0299 precluster_[i].ptOverNeighLocalMaxSum = (number > 1 ? 0.5 : 1.0) * rawet_[i];
0300 } break;
0301 }
0302 }
0303 }
0304
0305 clusterIndex_.fill(-1);
0306 clusters_.clear();
0307 unclustered_ = rawet_;
0308
0309 Cluster cluster;
0310 for (i = 0; i < ncells; ++i) {
0311 if (precluster_[i].ptLocalMax > 0) {
0312 float myet = rawet_[i];
0313 float tot = myet;
0314 float avg_eta = 0;
0315 float avg_phi = 0;
0316 cluster.clear();
0317 cluster.constituents.emplace_back(i, 1.0);
0318 for (const auto &ineigh : neighborCells_) {
0319 int ineighcell = grid_->neighbour(i, ineigh);
0320 if (ineighcell == -1)
0321 continue;
0322 float fracet = 0;
0323 switch (energyShareAlgo_) {
0324 case EnergyShareAlgo::Fractions:
0325 fracet = myet * precluster_.neigh(i, ineigh).ptOverNeighLocalMaxSum;
0326 break;
0327 case EnergyShareAlgo::None:
0328 fracet = precluster_.neigh(i, ineigh).ptOverNeighLocalMaxSum;
0329 break;
0330 case EnergyShareAlgo::Greedy:
0331 fracet = (myet == precluster_.neigh(i, ineigh).ptOverNeighLocalMaxSum ? rawet_.neigh(i, ineigh) : 0);
0332 break;
0333 case EnergyShareAlgo::Crude:
0334 fracet = precluster_.neigh(i, ineigh).ptOverNeighLocalMaxSum;
0335 break;
0336 }
0337 if (fracet == 0)
0338 continue;
0339 tot += fracet;
0340 cluster.constituents.emplace_back(ineighcell, fracet / rawet_.neigh(i, ineigh));
0341 if (energyWeightedPosition_) {
0342 avg_eta += fracet * (grid_->eta(ineighcell) - grid_->eta(i));
0343 avg_phi += fracet * deltaPhi(grid_->phi(ineighcell), grid_->phi(i));
0344 }
0345 }
0346 if (tot > minClusterEt_) {
0347 cluster.et = tot;
0348 unclustered_[i] = 0;
0349 for (const auto &ineigh : neighborCells_) {
0350 int ineighcell = grid_->neighbour(i, ineigh);
0351 if (ineighcell == -1)
0352 continue;
0353 unclustered_[ineighcell] = 0;
0354 }
0355 if (energyWeightedPosition_) {
0356 cluster.eta = grid_->eta(i) + avg_eta / tot;
0357 cluster.phi = grid_->phi(i) + avg_phi / tot;
0358
0359 cluster.phi = reco::reducePhiRange(cluster.phi);
0360 } else {
0361 cluster.eta = grid_->eta(i);
0362 cluster.phi = grid_->phi(i);
0363 }
0364 clusterIndex_[i] = clusters_.size();
0365 clusters_.push_back(cluster);
0366 }
0367 }
0368 }
0369 if (minEtToGrow_ > 0)
0370 grow();
0371 }
0372
0373 void l1tpf_calo::SingleCaloClusterer::grow() {
0374 int selneighs[4] = {1, 3, 4, 6};
0375 std::vector<int> toreset;
0376 for (Cluster &cluster : clusters_) {
0377 if (cluster.et > minEtToGrow_) {
0378 int i = cluster.constituents.front().first;
0379 for (int side = 0; side < 4; ++side) {
0380 int neigh = grid_->neighbour(i, selneighs[side]);
0381 if (neigh == -1)
0382 continue;
0383 for (int in = 0; in < 8; ++in) {
0384 int n2 = grid_->neighbour(neigh, in);
0385 if (n2 == -1)
0386 continue;
0387 cluster.et += unclustered_[n2];
0388 if (unclustered_[n2]) {
0389 cluster.constituents.emplace_back(n2, 1.0);
0390 toreset.push_back(n2);
0391 }
0392 }
0393 }
0394 }
0395 }
0396 for (int i : toreset)
0397 unclustered_[i] = 0;
0398 }
0399
0400 std::unique_ptr<l1t::PFClusterCollection> l1tpf_calo::SingleCaloClusterer::fetchCells(bool unclusteredOnly,
0401 float ptMin) const {
0402 auto ret = std::make_unique<l1t::PFClusterCollection>();
0403 const EtGrid &src = (unclusteredOnly ? unclustered_ : rawet_);
0404 const EtaPhiCenterGrid &eta_shift = eta_center_;
0405 const EtaPhiCenterGrid &phi_shift = phi_center_;
0406 l1tpf_calo::GridSelector selector = l1tpf_calo::GridSelector(etaBounds_, phiBounds_, maxClustersEtaPhi_);
0407 int totalClusters = 0;
0408 for (unsigned int i = 0, ncells = grid_->size(); i < ncells; ++i) {
0409 if (src[i] <= ptMin)
0410 continue;
0411 if ((unclusteredOnly == false) && (ptMin == 0)) {
0412 assert(cellKey_[i] == totalClusters);
0413 }
0414 totalClusters++;
0415 selector.fill(src[i], grid_->eta(i), grid_->phi(i), i);
0416 }
0417 std::vector<unsigned int> indices = selector.returnSorted();
0418 for (unsigned int ii = 0; ii < indices.size(); ii++) {
0419 unsigned int theIndex = indices[ii];
0420 ret->emplace_back(
0421 src[theIndex], grid_->eta(theIndex) + eta_shift[theIndex], grid_->phi(theIndex) + phi_shift[theIndex]);
0422 ret->back().setHwEta(grid_->ieta(theIndex));
0423 ret->back().setHwPhi(grid_->iphi(theIndex));
0424 }
0425 return ret;
0426 }
0427
0428 std::unique_ptr<l1t::PFClusterCollection> l1tpf_calo::SingleCaloClusterer::fetch(float ptMin) const {
0429 auto ret = std::make_unique<l1t::PFClusterCollection>();
0430 for (const Cluster &cluster : clusters_) {
0431 if (cluster.et > ptMin) {
0432 ret->emplace_back(cluster.et, cluster.eta, cluster.phi);
0433 }
0434 }
0435 return ret;
0436 }
0437
0438 std::unique_ptr<l1t::PFClusterCollection> l1tpf_calo::SingleCaloClusterer::fetch(
0439 const edm::OrphanHandle<l1t::PFClusterCollection> &cells, float ptMin) const {
0440 auto ret = std::make_unique<l1t::PFClusterCollection>();
0441 for (const Cluster &cluster : clusters_) {
0442 if (cluster.et > ptMin) {
0443 ret->emplace_back(cluster.et, cluster.eta, cluster.phi);
0444 for (const auto &pair : cluster.constituents) {
0445 edm::Ptr<l1t::PFCluster> ref(cells, cellKey_[pair.first]);
0446 ret->back().addConstituent(ref, pair.second);
0447 }
0448 }
0449 }
0450 return ret;
0451 }
0452
0453 l1tpf_calo::SimpleCaloLinkerBase::SimpleCaloLinkerBase(const edm::ParameterSet &pset,
0454 const SingleCaloClusterer &ecal,
0455 const SingleCaloClusterer &hcal)
0456 : grid_(getGrid(pset.getParameter<std::string>("grid"))),
0457 ecal_(ecal),
0458 hcal_(hcal),
0459 clusterIndex_(*grid_),
0460 clusters_(),
0461 etaBounds_(pset.getParameter<std::vector<double>>("etaBounds")),
0462 phiBounds_(pset.getParameter<std::vector<double>>("phiBounds")),
0463 maxClustersEtaPhi_(pset.getParameter<std::vector<unsigned int>>("maxClustersEtaPhi")),
0464 hoeCut_(pset.getParameter<double>("hoeCut")),
0465 minPhotonEt_(pset.getParameter<double>("minPhotonEt")),
0466 minHadronRawEt_(pset.getParameter<double>("minHadronRawEt")),
0467 minHadronEt_(pset.getParameter<double>("minHadronEt")),
0468 noEmInHGC_(pset.getParameter<bool>("noEmInHGC")) {
0469 if (grid_ != &ecal.raw().grid())
0470 throw cms::Exception("LogicError", "Inconsistent grid between ecal and linker\n");
0471 if (grid_ != &hcal.raw().grid())
0472 throw cms::Exception("LogicError", "Inconsistent grid between hcal and linker\n");
0473 if ((etaBounds_.size() - 1) * (phiBounds_.size() - 1) != maxClustersEtaPhi_.size()) {
0474 throw cms::Exception("Configuration")
0475 << "Size mismatch between eta/phi bounds and max clusters: " << (etaBounds_.size() - 1) << " x "
0476 << (phiBounds_.size() - 1) << " != " << maxClustersEtaPhi_.size() << "\n";
0477 }
0478 if (!std::is_sorted(etaBounds_.begin(), etaBounds_.end())) {
0479 throw cms::Exception("Configuration") << "etaBounds is not sorted\n";
0480 }
0481 if (!std::is_sorted(phiBounds_.begin(), phiBounds_.end())) {
0482 throw cms::Exception("Configuration") << "phiBounds is not sorted\n";
0483 }
0484 }
0485
0486 l1tpf_calo::SimpleCaloLinkerBase::~SimpleCaloLinkerBase() {}
0487
0488 std::unique_ptr<l1t::PFClusterCollection> l1tpf_calo::SimpleCaloLinkerBase::fetch() const {
0489 edm::OrphanHandle<l1t::PFClusterCollection> ecal, hcal;
0490 return fetch(ecal, hcal);
0491 }
0492
0493 std::unique_ptr<l1t::PFClusterCollection> l1tpf_calo::SimpleCaloLinkerBase::fetch(
0494 const edm::OrphanHandle<l1t::PFClusterCollection> &ecal,
0495 const edm::OrphanHandle<l1t::PFClusterCollection> &hcal) const {
0496 bool setRefs = (ecal.isValid() && hcal.isValid());
0497 auto ret = std::make_unique<l1t::PFClusterCollection>();
0498 l1tpf_calo::GridSelector selector = l1tpf_calo::GridSelector(etaBounds_, phiBounds_, maxClustersEtaPhi_);
0499 unsigned int index = 0;
0500 for (const CombinedCluster &cluster : clusters_) {
0501 index++;
0502 if (cluster.et > 0) {
0503 bool photon = (cluster.hcal_et < hoeCut_ * cluster.ecal_et);
0504 if (photon && noEmInHGC_) {
0505 if (std::abs(cluster.eta) > 1.5 && std::abs(cluster.eta) < 3.0) {
0506 continue;
0507 }
0508 }
0509 selector.fill(cluster.et, cluster.eta, cluster.phi, index - 1);
0510 }
0511 }
0512 std::vector<unsigned int> indices = selector.returnSorted();
0513 for (unsigned int ii = 0; ii < indices.size(); ii++) {
0514 unsigned int theIndex = indices[ii];
0515 const CombinedCluster &cluster = clusters_[theIndex];
0516 bool photon = (cluster.hcal_et < hoeCut_ * cluster.ecal_et);
0517 if (cluster.et > (photon ? minPhotonEt_ : minHadronEt_)) {
0518 ret->emplace_back(cluster.et,
0519 photon ? cluster.ecal_eta : cluster.eta,
0520 photon ? cluster.ecal_phi : cluster.phi,
0521 cluster.ecal_et > 0 ? std::max(cluster.et - cluster.ecal_et, 0.f) / cluster.ecal_et : -1,
0522 photon);
0523 if (setRefs) {
0524 for (const auto &pair : cluster.constituents) {
0525 assert(pair.first != 0);
0526 if (pair.first > 0) {
0527 ret->back().addConstituent(edm::Ptr<l1t::PFCluster>(hcal, +pair.first - 1), pair.second);
0528 } else {
0529 ret->back().addConstituent(edm::Ptr<l1t::PFCluster>(ecal, -pair.first + 1), pair.second);
0530 }
0531 }
0532 }
0533 }
0534 }
0535 return ret;
0536 }
0537
0538 l1tpf_calo::SimpleCaloLinker::SimpleCaloLinker(const edm::ParameterSet &pset,
0539 const SingleCaloClusterer &ecal,
0540 const SingleCaloClusterer &hcal)
0541 : SimpleCaloLinkerBase(pset, ecal, hcal), ecalToHCal_(*grid_) {}
0542
0543 l1tpf_calo::SimpleCaloLinker::~SimpleCaloLinker() {}
0544
0545 void l1tpf_calo::SimpleCaloLinker::clear() {
0546 clearBase();
0547 ecalToHCal_.clear();
0548 }
0549
0550 void l1tpf_calo::SimpleCaloLinker::run() {
0551 unsigned int i, ncells = grid_->size();
0552
0553 const EtGrid &hraw = hcal_.raw();
0554 const IndexGrid &ecals = ecal_.indexGrid();
0555 const IndexGrid &hcals = hcal_.indexGrid();
0556
0557
0558 ecalToHCal_.clear();
0559 for (i = 0; i < ncells; ++i) {
0560 if (ecals[i] >= 0) {
0561 if (hcals[i] >= 0) {
0562 ecalToHCal_[i].ptLocalMax = hcal_.cluster(i).et;
0563 } else {
0564 float tot = 0;
0565 for (int ineigh = 0; ineigh < 8; ++ineigh) {
0566 tot += hcal_.cluster(grid_->neighbour(i, ineigh)).et;
0567 }
0568 ecalToHCal_[i].ptOverNeighLocalMaxSum = tot ? ecal_.cluster(i).et / tot : 0;
0569 }
0570 }
0571 }
0572
0573 clusterIndex_.fill(-1);
0574 clusters_.clear();
0575 CombinedCluster cluster;
0576
0577 for (i = 0; i < ncells; ++i) {
0578 if (hcals[i] >= 0) {
0579 const Cluster &hcal = hcal_.cluster(i);
0580 cluster.clear();
0581 cluster.constituents.emplace_back(+i + 1, 1);
0582 if (ecalToHCal_[i].ptLocalMax > 0) {
0583
0584 const Cluster &ecal = ecal_.cluster(i);
0585 if (ecal.et + hcal.et > minHadronRawEt_) {
0586 cluster.ecal_et = ecal.et;
0587 cluster.hcal_et = hcal.et;
0588 cluster.et = cluster.ecal_et + cluster.hcal_et;
0589 float wecal = cluster.ecal_et / cluster.et, whcal = 1.0 - wecal;
0590 cluster.eta = ecal.eta * wecal + hcal.eta * whcal;
0591 cluster.phi = ecal.phi * wecal + hcal.phi * whcal;
0592 cluster.ecal_eta = cluster.eta;
0593 cluster.ecal_phi = cluster.phi;
0594
0595 cluster.phi = reco::reducePhiRange(cluster.phi);
0596 cluster.constituents.emplace_back(-i - 1, 1);
0597 }
0598 } else {
0599
0600 float myet = hcal.et;
0601 float etot = 0;
0602 float avg_eta = 0;
0603 float avg_phi = 0;
0604 for (int ineigh = 0; ineigh < 8; ++ineigh) {
0605 int ineighcell = grid_->neighbour(i, ineigh);
0606 if (ineighcell == -1)
0607 continue;
0608 float fracet = myet * ecalToHCal_.neigh(i, ineigh).ptOverNeighLocalMaxSum;
0609 if (fracet == 0)
0610 continue;
0611 etot += fracet;
0612 avg_eta += fracet * (grid_->eta(ineighcell) - grid_->eta(i));
0613 avg_phi += fracet * deltaPhi(grid_->phi(ineighcell), grid_->phi(i));
0614 cluster.constituents.emplace_back(-i - 1, fracet / ecal_.cluster(ineighcell).et);
0615 }
0616 if (myet + etot > minHadronRawEt_) {
0617 cluster.hcal_et = hcal.et;
0618 cluster.ecal_et = etot;
0619 cluster.et = myet + etot;
0620 cluster.eta = hcal.eta + avg_eta / cluster.et;
0621 cluster.phi = hcal.phi + avg_phi / cluster.et;
0622 cluster.ecal_eta = cluster.eta;
0623 cluster.ecal_phi = cluster.phi;
0624
0625 cluster.phi = reco::reducePhiRange(cluster.phi);
0626 }
0627 }
0628 if (cluster.et > 0) {
0629 clusterIndex_[i] = clusters_.size();
0630 clusters_.push_back(cluster);
0631 }
0632 }
0633 }
0634
0635
0636 for (i = 0; i < ncells; ++i) {
0637 if (ecals[i] >= 0 && ecalToHCal_[i].ptLocalMax == 0 && ecalToHCal_[i].ptOverNeighLocalMaxSum == 0) {
0638 cluster.clear();
0639 const Cluster &ecal = ecal_.cluster(i);
0640 cluster.ecal_et = ecal.et;
0641 cluster.hcal_et = hraw[i];
0642 cluster.et = cluster.ecal_et + cluster.hcal_et;
0643 cluster.eta = ecal.eta;
0644 cluster.phi = ecal.phi;
0645 cluster.constituents.emplace_back(-i - 1, 1);
0646 clusterIndex_[i] = clusters_.size();
0647 clusters_.push_back(cluster);
0648 }
0649 }
0650 }
0651
0652 l1tpf_calo::FlatCaloLinker::FlatCaloLinker(const edm::ParameterSet &pset,
0653 const SingleCaloClusterer &ecal,
0654 const SingleCaloClusterer &hcal)
0655 : SimpleCaloLinkerBase(pset, ecal, hcal), combClusterer_(pset) {}
0656
0657 l1tpf_calo::FlatCaloLinker::~FlatCaloLinker() {}
0658
0659 void l1tpf_calo::FlatCaloLinker::clear() {
0660 clearBase();
0661 combClusterer_.clear();
0662 }
0663
0664 void l1tpf_calo::FlatCaloLinker::run() {
0665 combClusterer_.clear();
0666
0667 const EtGrid &hraw = hcal_.raw();
0668 const EtGrid &eraw = ecal_.raw();
0669 combClusterer_.raw() = eraw;
0670 combClusterer_.raw() += hraw;
0671
0672 combClusterer_.run();
0673 clusterIndex_ = combClusterer_.indexGrid();
0674 const std::vector<Cluster> &clustersSrc = combClusterer_.clusters();
0675 unsigned int nclust = clustersSrc.size();
0676 clusters_.resize(nclust);
0677 for (unsigned int ic = 0; ic < nclust; ++ic) {
0678 const Cluster &src = clustersSrc[ic];
0679 CombinedCluster &dst = clusters_[ic];
0680 dst.et = src.et;
0681 dst.eta = src.eta;
0682 dst.phi = src.phi;
0683 dst.ecal_eta = src.eta;
0684 dst.ecal_phi = src.phi;
0685 dst.ecal_et = 0;
0686 dst.hcal_et = 0;
0687 for (const auto &pair : src.constituents) {
0688 if (eraw[pair.first]) {
0689 dst.ecal_et += pair.second * eraw[pair.first];
0690 dst.constituents.emplace_back(-pair.first - 1, pair.second);
0691 }
0692 if (hraw[pair.first]) {
0693 dst.hcal_et += pair.second * hraw[pair.first];
0694 dst.constituents.emplace_back(+pair.first + 1, pair.second);
0695 }
0696 }
0697 }
0698 }
0699
0700 l1tpf_calo::CombinedCaloLinker::CombinedCaloLinker(const edm::ParameterSet &pset,
0701 const SingleCaloClusterer &ecal,
0702 const SingleCaloClusterer &hcal)
0703 : SimpleCaloLinkerBase(pset, ecal, hcal), combClusterer_(pset) {}
0704
0705 l1tpf_calo::CombinedCaloLinker::~CombinedCaloLinker() {}
0706
0707 void l1tpf_calo::CombinedCaloLinker::clear() {
0708 clearBase();
0709 combClusterer_.clear();
0710 }
0711
0712 void l1tpf_calo::CombinedCaloLinker::run() {
0713 combClusterer_.clear();
0714
0715 const EtGrid &hraw = hcal_.raw();
0716 const EtGrid &eraw = ecal_.raw();
0717 const EtaPhiCenterGrid &eeta = ecal_.etaCenter();
0718 const EtaPhiCenterGrid &ephi = ecal_.phiCenter();
0719 combClusterer_.raw() = eraw;
0720 combClusterer_.raw() += hraw;
0721
0722 combClusterer_.run();
0723 clusterIndex_ = combClusterer_.indexGrid();
0724 const std::vector<Cluster> &clustersSrc = combClusterer_.clusters();
0725 unsigned int nclust = clustersSrc.size();
0726 clusters_.resize(nclust);
0727 for (unsigned int ic = 0; ic < nclust; ++ic) {
0728 const Cluster &src = clustersSrc[ic];
0729 CombinedCluster &dst = clusters_[ic];
0730 dst.et = src.et;
0731 dst.eta = src.eta;
0732 dst.phi = src.phi;
0733 dst.ecal_et = 0;
0734 dst.hcal_et = 0;
0735 float pt_max = 0.;
0736 float eta_ecal = 0.;
0737 float phi_ecal = 0.;
0738 for (const auto &pair : src.constituents) {
0739 if (eraw[pair.first]) {
0740 float ept = pair.second * eraw[pair.first];
0741 dst.ecal_et += ept;
0742 dst.constituents.emplace_back(-pair.first - 1, pair.second);
0743 if (ept > pt_max) {
0744 eta_ecal = eeta[pair.first];
0745 phi_ecal = ephi[pair.first];
0746 pt_max = ept;
0747 }
0748 }
0749 if (hraw[pair.first]) {
0750 dst.hcal_et += pair.second * hraw[pair.first];
0751 dst.constituents.emplace_back(+pair.first + 1, pair.second);
0752 }
0753 }
0754 dst.ecal_eta = eta_ecal;
0755 dst.ecal_phi = phi_ecal;
0756 }
0757 }
0758
0759 l1tpf_calo::GridSelector::GridSelector(std::vector<double> etaBounds,
0760 std::vector<double> phiBounds,
0761 std::vector<unsigned int> maxClusters)
0762 : etaBounds_(etaBounds),
0763 phiBounds_(phiBounds),
0764 maxClustersEtaPhi_(maxClusters),
0765 regionPtIndices_(!maxClusters.empty() ? maxClusters.size() : 1) {}
0766
0767 void l1tpf_calo::GridSelector::fill(float pt, float eta, float phi, unsigned int index) {
0768 if (!maxClustersEtaPhi_.empty()) {
0769 unsigned int etai = etaBounds_.size();
0770 for (unsigned int ie = 0; ie < etaBounds_.size() - 1; ie++) {
0771 if (eta >= etaBounds_[ie] && eta < etaBounds_[ie + 1]) {
0772 etai = ie;
0773 break;
0774 }
0775 }
0776 unsigned int phii = phiBounds_.size();
0777 for (unsigned int ip = 0; ip < phiBounds_.size() - 1; ip++) {
0778 if (phi >= phiBounds_[ip] && phi < phiBounds_[ip + 1]) {
0779 phii = ip;
0780 break;
0781 }
0782 }
0783 if (etai < etaBounds_.size() && phii < phiBounds_.size()) {
0784 regionPtIndices_[etai * (phiBounds_.size() - 1) + phii].emplace_back(pt, index);
0785 }
0786 } else {
0787 regionPtIndices_[0].emplace_back(pt, index);
0788 }
0789 }
0790
0791 std::vector<unsigned int> l1tpf_calo::GridSelector::returnSorted() {
0792 std::vector<unsigned int> indices;
0793 for (auto ®ionPtIndex : regionPtIndices_) {
0794 std::sort(regionPtIndex.begin(), regionPtIndex.end(), std::greater<std::pair<float, unsigned int>>());
0795 for (const auto &p : regionPtIndex) {
0796 indices.push_back(p.second);
0797 }
0798 }
0799 return indices;
0800 }
0801
0802 std::unique_ptr<l1tpf_calo::SimpleCaloLinkerBase> l1tpf_calo::makeCaloLinker(const edm::ParameterSet &pset,
0803 const SingleCaloClusterer &ecal,
0804 const SingleCaloClusterer &hcal) {
0805 const std::string &algo = pset.getParameter<std::string>("algo");
0806 if (algo == "simple") {
0807 return std::make_unique<l1tpf_calo::SimpleCaloLinker>(pset, ecal, hcal);
0808 } else if (algo == "flat") {
0809 return std::make_unique<l1tpf_calo::FlatCaloLinker>(pset, ecal, hcal);
0810 } else if (algo == "combined") {
0811 return std::make_unique<l1tpf_calo::CombinedCaloLinker>(pset, ecal, hcal);
0812 } else {
0813 throw cms::Exception("Configuration") << "Unsupported linker algo '" << algo << "'\n";
0814 }
0815 }