Back to home page

Project CMSSW displayed by LXR
 
 

    

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   // 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       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   // kill zeros. count non-zeros, for linking later
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   // pre-cluster step 1: at each cell, set the value equal to itself if it's a local maxima, zero otherwise
0217   // can be done in parallel on all cells
0218   for (i = 0; i < ncells; ++i) {
0219     if (rawet_[i] > seedEt_) {
0220       precluster_[i].ptLocalMax = rawet_[i];
0221       //// uncommment code below for debugging the clustering
0222       //printf("   candidate precluster pt %7.2f at %4d (ieta %+3d iphi %2d)\n",  rawet_[i], i, grid_->ieta(i), grid_->iphi(i));
0223       for (int ineigh = 0; ineigh <= 3; ++ineigh) {
0224         if (rawet_.neigh(i, ineigh) > rawet_[i])
0225           precluster_[i].ptLocalMax = 0;
0226         //// uncommment code below for debugging the clustering
0227         //int ncell = grid_->neighbour(i,ineigh);
0228         //if (ncell == -1) printf("   \t neigh %d is null\n", ineigh);
0229         //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);
0230       }
0231       for (int ineigh = 4; ineigh < 8; ++ineigh) {
0232         if (rawet_.neigh(i, ineigh) >= rawet_[i])
0233           precluster_[i].ptLocalMax = 0;
0234         //// uncommment code below for debugging the clustering
0235         //int ncell = grid_->neighbour(i,ineigh);
0236         //if (ncell == -1) printf("   \t neigh %d is null\n", ineigh);
0237         //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);
0238       }
0239     }
0240   }
0241   // pre-cluster step 2: compute information from neighbouring local max, for energy sharing purposes
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   // cluster: at each localMax cell, take itself plus the weighted contributions of the neighbours
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;  // skip dummy cells
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;  // skip dummy cells
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           // wrap around phi
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};  // -eta, -phi, +phi, +eta
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) {  // 1.5-3 = eta range of HGCal
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) {  // 1+hcal index
0459               ret->back().addConstituent(edm::Ptr<l1t::PFCluster>(hcal, +pair.first - 1), pair.second);
0460             } else {  // -1-ecal index
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   // for each ECal cluster, get the corresponding HCal cluster and the sum of the neighbour HCal clusters
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   // promote HCal clusters to final clusters
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         // direct linking is easy
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           // wrap around phi
0526           cluster.phi = reco::reduceRange(cluster.phi);
0527           cluster.constituents.emplace_back(-i - 1, 1);
0528         }
0529       } else {
0530         // sidewas linking is more annonying
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;  // skip dummy cells
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           // wrap around phi
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   // promote Unlinked ECal clusters to final clusters
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 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!