Phase2L1ParticleFlow/src/CaloClusterer.cc

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   // now link the cells
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   //// consistency check 1: check that find_cell works
0068   //// uncomment to check that there's no holes in the grid
0069   //for (float teta = 0; teta <= 5.0; teta += 0.02) {
0070   //    for (float tphi = -M_PI; tphi <= M_PI; tphi += 0.02) {
0071   //        find_cell(+teta, tphi);
0072   //        find_cell(-teta, tphi);
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;  // outside bounds
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);  // [-PI, PI]
0088   if (phi < 0)                   // then bring to [0, 2*PI]
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;  // fix corner case due to roundings etc
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   //// uncomment to check validity of derived coordinates
0100   //if (!valid_ieta_iphi(ieta,iphi)) {
0101   //    printf("Error in finding cell for eta %+7.4f phi %+7.4f, got ieta = %+3d iphi %2d which is not valid\n",
0102   //        eta, phi, ieta, iphi);
0103   //}
0104   assert(valid_ieta_iphi(ieta, iphi));
0105   int icell = ifind_cell(ieta, iphi);
0106   assert(icell != -1);
0107
0108   //// uncomment to check that the point is really in the cell
0109   //if (std::abs(eta - eta_[icell]) > 0.501*etaWidth_[icell] || std::abs(deltaPhi(phi, phi_[icell])) > 0.501*phiWidth_[icell]) {
0110   //    printf("Mismatch in finding cell for eta %+7.4f phi %+7.4f, got ieta = %+3d iphi %2d which has eta %+7.4f +- %.4f phi %+7.4f +- %.4f ; deta = %+7.4f dphi = %+7.4f\n",
0111   //        eta, phi, ieta, iphi, eta_[icell], etaWidth_[icell], phi_[icell], phiWidth_[icell], eta - eta_[icell], deltaPhi(phi, phi_[icell]));
0112   //}
0113   //assert(std::abs(eta - eta_[icell]) <= 0.5*etaWidth_[icell]);
0114   //assert(std::abs(deltaPhi(phi, phi_[icell])) <= 0.5*phiWidth_[icell]);
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");  //anything other than 3x3 is incompatible with grow() I think...
0200   } else {
0201     neighborCells_ = std::vector<int>({0, 1, 2, 3, 4, 5, 6, 7});  //default to 3x3
0202     //  in relative eta,phi: 5 = (+1, 0), 6 = (+1, 0), 7 = (+1,+1)
0203     //                       3 = ( 0,-1),              4 = ( 0,+1),
0204     //                       0 = (-1,-1), 1 = (-1, 0), 2 = (-1,+1),
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   // kill zeros. count non-zeros, for linking later
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   // pre-cluster step 1: at each cell, set the value equal to itself if it's a local maxima, zero otherwise
0245   // can be done in parallel on all cells
0246   for (i = 0; i < ncells; ++i) {
0247     if (rawet_[i] > seedEt_) {
0248       precluster_[i].ptLocalMax = rawet_[i];
0249       //// uncommment code below for debugging the clustering
0250       //printf("   candidate precluster pt %7.2f at %4d (ieta %+3d iphi %2d)\n",  rawet_[i], i, grid_->ieta(i), grid_->iphi(i));
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         //// uncommment code below for debugging the clustering
0257         //int ncell = grid_->neighbour(i,ineigh);
0258         //if (ncell == -1) printf("   \t neigh %d is null\n", ineigh);
0259         //else printf("   \t neigh %d at %4d (ieta %+3d iphi %2d) has pt %7.2f: comparison %1d \n", ineigh, ncell, grid_->ieta(ncell), grid_->iphi(ncell), rawet_[ncell], precluster_[i].ptLocalMax > 0);
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         //// uncommment code below for debugging the clustering
0267         //int ncell = grid_->neighbour(i,ineigh);
0268         //if (ncell == -1) printf("   \t neigh %d is null\n", ineigh);
0269         //else printf("   \t neigh %d at %4d (ieta %+3d iphi %2d) has pt %7.2f: comparison %1d \n", ineigh, ncell, grid_->ieta(ncell), grid_->iphi(ncell), rawet_[ncell], precluster_[i].ptLocalMax > 0);
0270       }
0271     }
0272   }
0273   // pre-cluster step 2: compute information from neighbouring local max, for energy sharing purposes
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   // cluster: at each localMax cell, take itself plus the weighted contributions of the neighbours
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;  // skip dummy cells
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;  // skip dummy cells
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           // wrap around phi
0359           cluster.phi = reco::reduceRange(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};  // -eta, -phi, +phi, +eta
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) {  // 1.5-3 = eta range of HGCal
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) {  // 1+hcal index
0527             ret->back().addConstituent(edm::Ptr<l1t::PFCluster>(hcal, +pair.first - 1), pair.second);
0528           } else {  // -1-ecal index
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   // for each ECal cluster, get the corresponding HCal cluster and the sum of the neighbour HCal clusters
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   // promote HCal clusters to final clusters
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         // direct linking is easy
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           // wrap around phi
0595           cluster.phi = reco::reduceRange(cluster.phi);
0596           cluster.constituents.emplace_back(-i - 1, 1);
0597         }
0598       } else {
0599         // sidewas linking is more annonying
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;  // skip dummy cells
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           // wrap around phi
0625           cluster.phi = reco::reduceRange(cluster.phi);
0626         }
0627       }
0628       if (cluster.et > 0) {
0629         clusterIndex_[i] = clusters_.size();
0630         clusters_.push_back(cluster);
0631       }
0632     }
0633   }
0634
0635   // promote Unlinked ECal clusters to final clusters
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 &regionPtIndex : 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 }