File indexing completed on 2021-02-14 13:30:27
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::reduceRange(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 precluster_(*grid_),
0171 clusterIndex_(*grid_),
0172 cellKey_(*grid_),
0173 clusters_(),
0174 nullCluster_(),
0175 zsEt_(pset.getParameter<double>("zsEt")),
0176 seedEt_(pset.getParameter<double>("seedEt")),
0177 minClusterEt_(pset.getParameter<double>("minClusterEt")),
0178 minEtToGrow_(pset.existsAs<double>("minEtToGrow") ? pset.getParameter<double>("minEtToGrow") : -1),
0179 energyWeightedPosition_(pset.getParameter<bool>("energyWeightedPosition")) {
0180 std::string energyShareAlgo = pset.getParameter<std::string>("energyShareAlgo");
0181 if (energyShareAlgo == "fractions")
0182 energyShareAlgo_ = EnergyShareAlgo::Fractions;
0183 else if (energyShareAlgo == "none")
0184 energyShareAlgo_ = EnergyShareAlgo::None;
0185 else if (energyShareAlgo == "greedy")
0186 energyShareAlgo_ = EnergyShareAlgo::Greedy;
0187 else if (energyShareAlgo == "crude")
0188 energyShareAlgo_ = EnergyShareAlgo::Crude;
0189 else
0190 throw cms::Exception("Configuration") << "Unsupported energyShareAlgo '" << energyShareAlgo << "'\n";
0191 }
0192
0193 l1tpf_calo::SingleCaloClusterer::~SingleCaloClusterer() {}
0194
0195 void l1tpf_calo::SingleCaloClusterer::clear() {
0196 rawet_.zero();
0197 clusters_.clear();
0198 clusterIndex_.fill(-1);
0199 }
0200
0201 void l1tpf_calo::SingleCaloClusterer::run() {
0202 unsigned int i, ncells = grid_->size();
0203
0204
0205 cellKey_.fill(-1);
0206 int key = 0;
0207 for (i = 0; i < ncells; ++i) {
0208 if (rawet_[i] < zsEt_) {
0209 rawet_[i] = 0;
0210 } else {
0211 cellKey_[i] = key++;
0212 }
0213 }
0214
0215 precluster_.clear();
0216
0217
0218 for (i = 0; i < ncells; ++i) {
0219 if (rawet_[i] > seedEt_) {
0220 precluster_[i].ptLocalMax = rawet_[i];
0221
0222
0223 for (int ineigh = 0; ineigh <= 3; ++ineigh) {
0224 if (rawet_.neigh(i, ineigh) > rawet_[i])
0225 precluster_[i].ptLocalMax = 0;
0226
0227
0228
0229
0230 }
0231 for (int ineigh = 4; ineigh < 8; ++ineigh) {
0232 if (rawet_.neigh(i, ineigh) >= rawet_[i])
0233 precluster_[i].ptLocalMax = 0;
0234
0235
0236
0237
0238 }
0239 }
0240 }
0241
0242 for (i = 0; i < ncells; ++i) {
0243 if (precluster_[i].ptLocalMax == 0) {
0244 switch (energyShareAlgo_) {
0245 case EnergyShareAlgo::Fractions: {
0246 float tot = 0;
0247 for (int ineigh = 0; ineigh < 8; ++ineigh) {
0248 tot += precluster_.neigh(i, ineigh).ptLocalMax;
0249 }
0250 precluster_[i].ptOverNeighLocalMaxSum = tot ? rawet_[i] / tot : 0;
0251 } break;
0252 case EnergyShareAlgo::None:
0253 precluster_[i].ptOverNeighLocalMaxSum = rawet_[i];
0254 break;
0255 case EnergyShareAlgo::Greedy: {
0256 float maxet = 0;
0257 for (int ineigh = 0; ineigh < 8; ++ineigh) {
0258 maxet = std::max(maxet, precluster_.neigh(i, ineigh).ptLocalMax);
0259 }
0260 precluster_[i].ptOverNeighLocalMaxSum = maxet;
0261 } break;
0262 case EnergyShareAlgo::Crude: {
0263 int number = 0;
0264 for (int ineigh = 0; ineigh < 8; ++ineigh) {
0265 number += (precluster_.neigh(i, ineigh).ptLocalMax > 0);
0266 }
0267 precluster_[i].ptOverNeighLocalMaxSum = (number > 1 ? 0.5 : 1.0) * rawet_[i];
0268 } break;
0269 }
0270 }
0271 }
0272
0273 clusterIndex_.fill(-1);
0274 clusters_.clear();
0275 unclustered_ = rawet_;
0276
0277 Cluster cluster;
0278 for (i = 0; i < ncells; ++i) {
0279 if (precluster_[i].ptLocalMax > 0) {
0280 float myet = rawet_[i];
0281 float tot = myet;
0282 float avg_eta = 0;
0283 float avg_phi = 0;
0284 cluster.clear();
0285 cluster.constituents.emplace_back(i, 1.0);
0286 for (int ineigh = 0; ineigh < 8; ++ineigh) {
0287 int ineighcell = grid_->neighbour(i, ineigh);
0288 if (ineighcell == -1)
0289 continue;
0290 float fracet = 0;
0291 switch (energyShareAlgo_) {
0292 case EnergyShareAlgo::Fractions:
0293 fracet = myet * precluster_.neigh(i, ineigh).ptOverNeighLocalMaxSum;
0294 break;
0295 case EnergyShareAlgo::None:
0296 fracet = precluster_.neigh(i, ineigh).ptOverNeighLocalMaxSum;
0297 break;
0298 case EnergyShareAlgo::Greedy:
0299 fracet = (myet == precluster_.neigh(i, ineigh).ptOverNeighLocalMaxSum ? rawet_.neigh(i, ineigh) : 0);
0300 break;
0301 case EnergyShareAlgo::Crude:
0302 fracet = precluster_.neigh(i, ineigh).ptOverNeighLocalMaxSum;
0303 break;
0304 }
0305 if (fracet == 0)
0306 continue;
0307 tot += fracet;
0308 cluster.constituents.emplace_back(ineighcell, fracet / rawet_.neigh(i, ineigh));
0309 if (energyWeightedPosition_) {
0310 avg_eta += fracet * (grid_->eta(ineighcell) - grid_->eta(i));
0311 avg_phi += fracet * deltaPhi(grid_->phi(ineighcell), grid_->phi(i));
0312 }
0313 }
0314 if (tot > minClusterEt_) {
0315 cluster.et = tot;
0316 unclustered_[i] = 0;
0317 for (int ineigh = 0; ineigh < 8; ++ineigh) {
0318 int ineighcell = grid_->neighbour(i, ineigh);
0319 if (ineighcell == -1)
0320 continue;
0321 unclustered_[ineighcell] = 0;
0322 }
0323 if (energyWeightedPosition_) {
0324 cluster.eta = grid_->eta(i) + avg_eta / tot;
0325 cluster.phi = grid_->phi(i) + avg_phi / tot;
0326
0327 cluster.phi = reco::reduceRange(cluster.phi);
0328 } else {
0329 cluster.eta = grid_->eta(i);
0330 cluster.phi = grid_->phi(i);
0331 }
0332 clusterIndex_[i] = clusters_.size();
0333 clusters_.push_back(cluster);
0334 }
0335 }
0336 }
0337 if (minEtToGrow_ > 0)
0338 grow();
0339 }
0340
0341 void l1tpf_calo::SingleCaloClusterer::grow() {
0342 int selneighs[4] = {1, 3, 4, 6};
0343 std::vector<int> toreset;
0344 for (Cluster &cluster : clusters_) {
0345 if (cluster.et > minEtToGrow_) {
0346 int i = cluster.constituents.front().first;
0347 for (int side = 0; side < 4; ++side) {
0348 int neigh = grid_->neighbour(i, selneighs[side]);
0349 if (neigh == -1)
0350 continue;
0351 for (int in = 0; in < 8; ++in) {
0352 int n2 = grid_->neighbour(neigh, in);
0353 if (n2 == -1)
0354 continue;
0355 cluster.et += unclustered_[n2];
0356 if (unclustered_[n2]) {
0357 cluster.constituents.emplace_back(n2, 1.0);
0358 toreset.push_back(n2);
0359 }
0360 }
0361 }
0362 }
0363 }
0364 for (int i : toreset)
0365 unclustered_[i] = 0;
0366 }
0367
0368 std::unique_ptr<l1t::PFClusterCollection> l1tpf_calo::SingleCaloClusterer::fetchCells(bool unclusteredOnly,
0369 float ptMin) const {
0370 auto ret = std::make_unique<l1t::PFClusterCollection>();
0371 const EtGrid &src = (unclusteredOnly ? unclustered_ : rawet_);
0372 for (unsigned int i = 0, ncells = grid_->size(); i < ncells; ++i) {
0373 if (src[i] <= ptMin)
0374 continue;
0375 if ((unclusteredOnly == false) && (ptMin == 0)) {
0376 assert(cellKey_[i] == int(ret->size()));
0377 }
0378 ret->emplace_back(src[i], grid_->eta(i), grid_->phi(i));
0379 ret->back().setHwEta(grid_->ieta(i));
0380 ret->back().setHwPhi(grid_->iphi(i));
0381 }
0382 return ret;
0383 }
0384
0385 std::unique_ptr<l1t::PFClusterCollection> l1tpf_calo::SingleCaloClusterer::fetch(float ptMin) const {
0386 auto ret = std::make_unique<l1t::PFClusterCollection>();
0387 for (const Cluster &cluster : clusters_) {
0388 if (cluster.et > ptMin) {
0389 ret->emplace_back(cluster.et, cluster.eta, cluster.phi);
0390 }
0391 }
0392 return ret;
0393 }
0394
0395 std::unique_ptr<l1t::PFClusterCollection> l1tpf_calo::SingleCaloClusterer::fetch(
0396 const edm::OrphanHandle<l1t::PFClusterCollection> &cells, float ptMin) const {
0397 auto ret = std::make_unique<l1t::PFClusterCollection>();
0398 for (const Cluster &cluster : clusters_) {
0399 if (cluster.et > ptMin) {
0400 ret->emplace_back(cluster.et, cluster.eta, cluster.phi);
0401 for (const auto &pair : cluster.constituents) {
0402 edm::Ptr<l1t::PFCluster> ref(cells, cellKey_[pair.first]);
0403 ret->back().addConstituent(ref, pair.second);
0404 }
0405 }
0406 }
0407 return ret;
0408 }
0409
0410 l1tpf_calo::SimpleCaloLinkerBase::SimpleCaloLinkerBase(const edm::ParameterSet &pset,
0411 const SingleCaloClusterer &ecal,
0412 const SingleCaloClusterer &hcal)
0413 : grid_(getGrid(pset.getParameter<std::string>("grid"))),
0414 ecal_(ecal),
0415 hcal_(hcal),
0416 clusterIndex_(*grid_),
0417 clusters_(),
0418 hoeCut_(pset.getParameter<double>("hoeCut")),
0419 minPhotonEt_(pset.getParameter<double>("minPhotonEt")),
0420 minHadronRawEt_(pset.getParameter<double>("minHadronRawEt")),
0421 minHadronEt_(pset.getParameter<double>("minHadronEt")),
0422 noEmInHGC_(pset.getParameter<bool>("noEmInHGC")) {
0423 if (grid_ != &ecal.raw().grid())
0424 throw cms::Exception("LogicError", "Inconsistent grid between ecal and linker\n");
0425 if (grid_ != &hcal.raw().grid())
0426 throw cms::Exception("LogicError", "Inconsistent grid between hcal and linker\n");
0427 }
0428
0429 l1tpf_calo::SimpleCaloLinkerBase::~SimpleCaloLinkerBase() {}
0430
0431 std::unique_ptr<l1t::PFClusterCollection> l1tpf_calo::SimpleCaloLinkerBase::fetch() const {
0432 edm::OrphanHandle<l1t::PFClusterCollection> ecal, hcal;
0433 return fetch(ecal, hcal);
0434 }
0435
0436 std::unique_ptr<l1t::PFClusterCollection> l1tpf_calo::SimpleCaloLinkerBase::fetch(
0437 const edm::OrphanHandle<l1t::PFClusterCollection> &ecal,
0438 const edm::OrphanHandle<l1t::PFClusterCollection> &hcal) const {
0439 bool setRefs = (ecal.isValid() && hcal.isValid());
0440 auto ret = std::make_unique<l1t::PFClusterCollection>();
0441 for (const CombinedCluster &cluster : clusters_) {
0442 if (cluster.et > 0) {
0443 bool photon = (cluster.hcal_et < hoeCut_ * cluster.ecal_et);
0444 if (photon && noEmInHGC_) {
0445 if (std::abs(cluster.eta) > 1.5 && std::abs(cluster.eta) < 3.0) {
0446 continue;
0447 }
0448 }
0449 if (cluster.et > (photon ? minPhotonEt_ : minHadronEt_)) {
0450 ret->emplace_back(cluster.et,
0451 cluster.eta,
0452 cluster.phi,
0453 cluster.ecal_et > 0 ? std::max(cluster.et - cluster.ecal_et, 0.f) / cluster.ecal_et : -1,
0454 photon);
0455 if (setRefs) {
0456 for (const auto &pair : cluster.constituents) {
0457 assert(pair.first != 0);
0458 if (pair.first > 0) {
0459 ret->back().addConstituent(edm::Ptr<l1t::PFCluster>(hcal, +pair.first - 1), pair.second);
0460 } else {
0461 ret->back().addConstituent(edm::Ptr<l1t::PFCluster>(ecal, -pair.first + 1), pair.second);
0462 }
0463 }
0464 }
0465 }
0466 }
0467 }
0468 return ret;
0469 }
0470
0471 l1tpf_calo::SimpleCaloLinker::SimpleCaloLinker(const edm::ParameterSet &pset,
0472 const SingleCaloClusterer &ecal,
0473 const SingleCaloClusterer &hcal)
0474 : SimpleCaloLinkerBase(pset, ecal, hcal), ecalToHCal_(*grid_) {}
0475
0476 l1tpf_calo::SimpleCaloLinker::~SimpleCaloLinker() {}
0477
0478 void l1tpf_calo::SimpleCaloLinker::clear() {
0479 clearBase();
0480 ecalToHCal_.clear();
0481 }
0482
0483 void l1tpf_calo::SimpleCaloLinker::run() {
0484 unsigned int i, ncells = grid_->size();
0485
0486 const EtGrid &hraw = hcal_.raw();
0487 const IndexGrid &ecals = ecal_.indexGrid();
0488 const IndexGrid &hcals = hcal_.indexGrid();
0489
0490
0491 ecalToHCal_.clear();
0492 for (i = 0; i < ncells; ++i) {
0493 if (ecals[i] >= 0) {
0494 if (hcals[i] >= 0) {
0495 ecalToHCal_[i].ptLocalMax = hcal_.cluster(i).et;
0496 } else {
0497 float tot = 0;
0498 for (int ineigh = 0; ineigh < 8; ++ineigh) {
0499 tot += hcal_.cluster(grid_->neighbour(i, ineigh)).et;
0500 }
0501 ecalToHCal_[i].ptOverNeighLocalMaxSum = tot ? ecal_.cluster(i).et / tot : 0;
0502 }
0503 }
0504 }
0505
0506 clusterIndex_.fill(-1);
0507 clusters_.clear();
0508 CombinedCluster cluster;
0509
0510 for (i = 0; i < ncells; ++i) {
0511 if (hcals[i] >= 0) {
0512 const Cluster &hcal = hcal_.cluster(i);
0513 cluster.clear();
0514 cluster.constituents.emplace_back(+i + 1, 1);
0515 if (ecalToHCal_[i].ptLocalMax > 0) {
0516
0517 const Cluster &ecal = ecal_.cluster(i);
0518 if (ecal.et + hcal.et > minHadronRawEt_) {
0519 cluster.ecal_et = ecal.et;
0520 cluster.hcal_et = hcal.et;
0521 cluster.et = cluster.ecal_et + cluster.hcal_et;
0522 float wecal = cluster.ecal_et / cluster.et, whcal = 1.0 - wecal;
0523 cluster.eta = ecal.eta * wecal + hcal.eta * whcal;
0524 cluster.phi = ecal.phi * wecal + hcal.phi * whcal;
0525
0526 cluster.phi = reco::reduceRange(cluster.phi);
0527 cluster.constituents.emplace_back(-i - 1, 1);
0528 }
0529 } else {
0530
0531 float myet = hcal.et;
0532 float etot = 0;
0533 float avg_eta = 0;
0534 float avg_phi = 0;
0535 for (int ineigh = 0; ineigh < 8; ++ineigh) {
0536 int ineighcell = grid_->neighbour(i, ineigh);
0537 if (ineighcell == -1)
0538 continue;
0539 float fracet = myet * ecalToHCal_.neigh(i, ineigh).ptOverNeighLocalMaxSum;
0540 if (fracet == 0)
0541 continue;
0542 etot += fracet;
0543 avg_eta += fracet * (grid_->eta(ineighcell) - grid_->eta(i));
0544 avg_phi += fracet * deltaPhi(grid_->phi(ineighcell), grid_->phi(i));
0545 cluster.constituents.emplace_back(-i - 1, fracet / ecal_.cluster(ineighcell).et);
0546 }
0547 if (myet + etot > minHadronRawEt_) {
0548 cluster.hcal_et = hcal.et;
0549 cluster.ecal_et = etot;
0550 cluster.et = myet + etot;
0551 cluster.eta = hcal.eta + avg_eta / cluster.et;
0552 cluster.phi = hcal.phi + avg_phi / cluster.et;
0553
0554 cluster.phi = reco::reduceRange(cluster.phi);
0555 }
0556 }
0557 if (cluster.et > 0) {
0558 clusterIndex_[i] = clusters_.size();
0559 clusters_.push_back(cluster);
0560 }
0561 }
0562 }
0563
0564
0565 for (i = 0; i < ncells; ++i) {
0566 if (ecals[i] >= 0 && ecalToHCal_[i].ptLocalMax == 0 && ecalToHCal_[i].ptOverNeighLocalMaxSum == 0) {
0567 cluster.clear();
0568 const Cluster &ecal = ecal_.cluster(i);
0569 cluster.ecal_et = ecal.et;
0570 cluster.hcal_et = hraw[i];
0571 cluster.et = cluster.ecal_et + cluster.hcal_et;
0572 cluster.eta = ecal.eta;
0573 cluster.phi = ecal.phi;
0574 cluster.constituents.emplace_back(-i - 1, 1);
0575 clusterIndex_[i] = clusters_.size();
0576 clusters_.push_back(cluster);
0577 }
0578 }
0579 }
0580
0581 l1tpf_calo::FlatCaloLinker::FlatCaloLinker(const edm::ParameterSet &pset,
0582 const SingleCaloClusterer &ecal,
0583 const SingleCaloClusterer &hcal)
0584 : SimpleCaloLinkerBase(pset, ecal, hcal), combClusterer_(pset) {}
0585
0586 l1tpf_calo::FlatCaloLinker::~FlatCaloLinker() {}
0587
0588 void l1tpf_calo::FlatCaloLinker::clear() {
0589 clearBase();
0590 combClusterer_.clear();
0591 }
0592
0593 void l1tpf_calo::FlatCaloLinker::run() {
0594 combClusterer_.clear();
0595
0596 const EtGrid &hraw = hcal_.raw();
0597 const EtGrid &eraw = ecal_.raw();
0598 combClusterer_.raw() = eraw;
0599 combClusterer_.raw() += hraw;
0600
0601 combClusterer_.run();
0602 clusterIndex_ = combClusterer_.indexGrid();
0603 const std::vector<Cluster> &clustersSrc = combClusterer_.clusters();
0604 unsigned int nclust = clustersSrc.size();
0605 clusters_.resize(nclust);
0606 for (unsigned int ic = 0; ic < nclust; ++ic) {
0607 const Cluster &src = clustersSrc[ic];
0608 CombinedCluster &dst = clusters_[ic];
0609 dst.et = src.et;
0610 dst.eta = src.eta;
0611 dst.phi = src.phi;
0612 dst.ecal_et = 0;
0613 dst.hcal_et = 0;
0614 for (const auto &pair : src.constituents) {
0615 if (eraw[pair.first]) {
0616 dst.ecal_et += pair.second * eraw[pair.first];
0617 dst.constituents.emplace_back(-pair.first - 1, pair.second);
0618 }
0619 if (hraw[pair.first]) {
0620 dst.hcal_et += pair.second * hraw[pair.first];
0621 dst.constituents.emplace_back(+pair.first + 1, pair.second);
0622 }
0623 }
0624 }
0625 }
0626
0627 std::unique_ptr<l1tpf_calo::SimpleCaloLinkerBase> l1tpf_calo::makeCaloLinker(const edm::ParameterSet &pset,
0628 const SingleCaloClusterer &ecal,
0629 const SingleCaloClusterer &hcal) {
0630 const std::string &algo = pset.getParameter<std::string>("algo");
0631 if (algo == "simple") {
0632 return std::make_unique<l1tpf_calo::SimpleCaloLinker>(pset, ecal, hcal);
0633 } else if (algo == "flat") {
0634 return std::make_unique<l1tpf_calo::FlatCaloLinker>(pset, ecal, hcal);
0635 } else {
0636 throw cms::Exception("Configuration") << "Unsupported linker algo '" << algo << "'\n";
0637 }
0638 }