PFClusterProducer/plugins/RBXAndHPDCleaner.cc

0001 #include "DataFormats/HcalDetId/interface/HcalDetId.h"
0002 #include "FWCore/Framework/interface/ConsumesCollector.h"
0003 #include "RecoParticleFlow/PFClusterProducer/interface/RecHitTopologicalCleanerBase.h"
0004
0005 #include <cmath>
0006 #include <unordered_map>
0007
0008 class RBXAndHPDCleaner : public RecHitTopologicalCleanerBase {
0009 public:
0010   RBXAndHPDCleaner(const edm::ParameterSet& conf, edm::ConsumesCollector& cc)
0011       : RecHitTopologicalCleanerBase(conf, cc) {}
0012   RBXAndHPDCleaner(const RBXAndHPDCleaner&) = delete;
0013   RBXAndHPDCleaner& operator=(const RBXAndHPDCleaner&) = delete;
0014
0015   void clean(const edm::Handle<reco::PFRecHitCollection>& input, std::vector<bool>& mask) override;
0016
0017 private:
0018   std::unordered_map<int, std::vector<unsigned> > _hpds, _rbxs;
0019 };
0020
0021 DEFINE_EDM_PLUGIN(RecHitTopologicalCleanerFactory, RBXAndHPDCleaner, "RBXAndHPDCleaner");
0022
0023 namespace {
0024   bool greaterByEnergy(const std::pair<unsigned, double>& a, const std::pair<unsigned, double>& b) {
0025     return a.second > b.second;
0026   }
0027 }  // namespace
0028
0029 void RBXAndHPDCleaner::clean(const edm::Handle<reco::PFRecHitCollection>& input, std::vector<bool>& mask) {
0030   _hpds.clear();
0031   _rbxs.clear();
0032   //need to run over energy sorted rechits
0033   std::vector<std::pair<unsigned, double> > ordered_hits;
0034   ordered_hits.reserve(input->size());
0035   for (unsigned i = 0; i < input->size(); ++i) {
0036     std::pair<unsigned, double> val = std::make_pair(i, input->at(i).energy());
0037     auto pos = std::upper_bound(ordered_hits.begin(), ordered_hits.end(), val, greaterByEnergy);
0038     ordered_hits.insert(pos, val);
0039   }
0040
0041   for (const auto& idx_e : ordered_hits) {
0042     if (!mask[idx_e.first])
0043       continue;
0044     const unsigned idx = idx_e.first;
0045     const reco::PFRecHit& rechit = input->at(idx);
0046     int layer = rechit.layer();
0047     if (layer != PFLayer::HCAL_BARREL1 && layer != PFLayer::HCAL_ENDCAP) {
0048       break;
0049     }
0050     HcalDetId theHcalDetId(rechit.detId());
0051     int ieta = theHcalDetId.ieta();
0052     int iphi = theHcalDetId.iphi();
0053     int ihpd = 0, irbx = 0;
0054     switch (layer) {
0055       case PFLayer::HCAL_BARREL1:
0056         ihpd = (ieta < 0 ? -iphi : iphi);
0057         irbx = (ieta < 0 ? -(iphi + 5) / 4 : (iphi + 5) / 4);
0058         break;
0059       case PFLayer::HCAL_ENDCAP:
0060         ihpd = (ieta < 0 ? -(iphi + 1) / 2 - 100 : (iphi + 1) / 2 + 100);
0061         irbx = (ieta < 0 ? -(iphi + 5) / 4 - 20 : (iphi + 5) / 4 + 20);
0062         break;
0063       default:
0064         break;
0065     }
0066     switch (std::abs(irbx)) {
0067       case 19:
0068         irbx = (irbx < 0 ? -1 : 1);
0069         break;
0070       case 39:
0071         irbx = (irbx < 0 ? -21 : 21);
0072         break;
0073       default:
0074         break;
0075     }
0076     _hpds[ihpd].push_back(idx);
0077     _rbxs[irbx].push_back(idx);
0078   }
0079   // loop on the rbx's we found and clean RBX's with tons of rechits
0080   // and lots of energy
0081   std::unordered_map<int, std::vector<unsigned> > theHPDs;
0082   std::unordered_multimap<double, unsigned> theEnergies;
0083   for (const auto& itrbx : _rbxs) {
0084     if ((std::abs(itrbx.first) < 20 && itrbx.second.size() > 30) ||
0085         (std::abs(itrbx.first) > 20 && itrbx.second.size() > 30)) {
0086       const std::vector<unsigned>& rechits = itrbx.second;
0087       theHPDs.clear();
0088       theEnergies.clear();
0089       int nSeeds0 = rechits.size();
0090       for (unsigned jh = 0; jh < rechits.size(); ++jh) {
0091         const reco::PFRecHit& rechit = (*input)[jh];
0092         // check if rechit is a seed
0093         unsigned nN = 0;  // neighbours over threshold
0094         bool isASeed = true;
0095         auto const& neighbours4 = rechit.neighbours4();
0096         for (auto k : neighbours4) {
0097           auto const& neighbour = (*input)[k];
0098           if (neighbour.energy() > rechit.energy()) {
0099             --nSeeds0;
0100             isASeed = false;
0101             break;
0102           } else {
0103             if (neighbour.energy() > 0.4)
0104               ++nN;
0105           }
0106         }
0107         if (isASeed && !nN)
0108           --nSeeds0;
0109
0110         HcalDetId theHcalDetId(rechit.detId());
0111         int iphi = theHcalDetId.iphi();
0112         switch (rechit.layer()) {
0113           case PFLayer::HCAL_BARREL1:
0114             theHPDs[iphi].push_back(rechits[jh]);
0115             break;
0116           case PFLayer::HCAL_ENDCAP:
0117             theHPDs[(iphi - 1) / 2].push_back(rechits[jh]);
0118             break;
0119           default:
0120             break;
0121         }
0122         const double rhenergy = rechit.energy();
0123         theEnergies.emplace(rhenergy, rechits[jh]);
0124       }
0125       if (nSeeds0 > 6) {
0126         unsigned nHPD15 = 0;
0127         for (const auto& itHPD : theHPDs) {
0128           int hpdN = itHPD.first;
0129           const std::vector<unsigned>& hpdHits = itHPD.second;
0130           if ((std::abs(hpdN) < 100 && hpdHits.size() > 14) || (std::abs(hpdN) > 100 && hpdHits.size() > 14))
0131             ++nHPD15;
0132         }
0133         if (nHPD15 > 1) {
0134           unsigned nn = 0;
0135           double threshold = 1.0;
0136           for (const auto& itEn : theEnergies) {
0137             ++nn;
0138             if (nn < 5) {
0139               mask[itEn.second] = false;
0140             } else if (nn == 5) {
0141               threshold = itEn.first * 5;
0142               mask[itEn.second] = false;
0143             } else {
0144               if (itEn.first < threshold)
0145                 mask[itEn.second] = false;
0146             }
0147           }
0148         }
0149       }
0150     }
0151   }
0152
0153   // loop on the HPDs
0154   std::unordered_map<int, std::vector<unsigned> >::iterator neighbour1;
0155   std::unordered_map<int, std::vector<unsigned> >::iterator neighbour2;
0156   std::unordered_map<int, std::vector<unsigned> >::iterator neighbour0;
0157   std::unordered_map<int, std::vector<unsigned> >::iterator neighbour3;
0158   unsigned size1 = 0, size2 = 0;
0159   for (const auto& ithpd : _hpds) {
0160     const std::vector<unsigned>& rechits = ithpd.second;
0161     theEnergies.clear();
0162     for (const unsigned rhidx : rechits) {
0163       const reco::PFRecHit& rechit = input->at(rhidx);
0164       theEnergies.emplace(rechit.energy(), rhidx);
0165     }
0166
0167     const int thehpd = ithpd.first;
0168     switch (std::abs(thehpd)) {
0169       case 1:
0170         neighbour1 = (thehpd > 0 ? _hpds.find(72) : _hpds.find(-72));
0171         break;
0172       case 72:
0173         neighbour2 = (thehpd > 0 ? _hpds.find(1) : _hpds.find(-1));
0174         break;
0175       case 101:
0176         neighbour1 = (thehpd > 0 ? _hpds.find(136) : _hpds.find(-136));
0177         break;
0178       case 136:
0179         neighbour2 = (thehpd > 0 ? _hpds.find(101) : _hpds.find(-101));
0180         break;
0181       default:
0182         neighbour1 = (thehpd > 0 ? _hpds.find(thehpd - 1) : _hpds.find(thehpd + 1));
0183         neighbour2 = (thehpd > 0 ? _hpds.find(thehpd + 1) : _hpds.find(thehpd - 1));
0184         break;
0185     }
0186     if (neighbour1 != _hpds.end()) {
0187       const int nb1 = neighbour1->first;
0188       switch (std::abs(nb1)) {
0189         case 1:
0190           neighbour0 = (nb1 > 0 ? _hpds.find(72) : _hpds.find(-72));
0191           break;
0192         case 101:
0193           neighbour0 = (nb1 > 0 ? _hpds.find(136) : _hpds.find(-136));
0194           break;
0195         default:
0196           neighbour0 = (nb1 > 0 ? _hpds.find(nb1 - 1) : _hpds.find(nb1 + 1));
0197           break;
0198       }
0199     } else {
0200       neighbour0 = _hpds.end();
0201     }
0202
0203     if (neighbour2 != _hpds.end()) {
0204       const int nb2 = neighbour2->first;
0205       switch (std::abs(nb2)) {
0206         case 72:
0207           neighbour3 = (nb2 > 0 ? _hpds.find(1) : _hpds.find(-1));
0208           break;
0209         case 136:
0210           neighbour3 = (nb2 > 0 ? _hpds.find(101) : _hpds.find(-101));
0211           break;
0212         default:
0213           neighbour3 = (nb2 > 0 ? _hpds.find(nb2 + 1) : _hpds.find(nb2 - 1));
0214           break;
0215       }
0216     } else {
0217       neighbour3 = _hpds.end();
0218     }
0219
0220     size1 = neighbour1 != _hpds.end() ? neighbour1->second.size() : 0;
0221     size2 = neighbour2 != _hpds.end() ? neighbour2->second.size() : 0;
0222     if (size1 > 10) {
0223       if ((abs(neighbour1->first) > 100 && neighbour1->second.size() > 15) ||
0224           (abs(neighbour1->first) < 100 && neighbour1->second.size() > 12))
0225         size1 = neighbour0 != _hpds.end() ? neighbour0->second.size() : 0;
0226     }
0227     if (size2 > 10) {
0228       if ((abs(neighbour2->first) > 100 && neighbour2->second.size() > 15) ||
0229           (abs(neighbour2->first) < 100 && neighbour2->second.size() > 12))
0230         size2 = neighbour3 != _hpds.end() ? neighbour3->second.size() : 0;
0231     }
0232     if ((std::abs(ithpd.first) > 100 && ithpd.second.size() > 15) ||
0233         (std::abs(ithpd.first) < 100 && ithpd.second.size() > 12)) {
0234       if ((double)(size1 + size2) / (float)ithpd.second.size() < 1.0) {
0235         unsigned nn = 0;
0236         double threshold = 1.0;
0237         for (const auto& itEn : theEnergies) {
0238           if (nn < 5) {
0239             mask[itEn.second] = false;
0240           } else if (nn == 5) {
0241             threshold = itEn.first * 2.5;
0242             mask[itEn.second] = false;
0243           } else {
0244             if (itEn.first < threshold)
0245               mask[itEn.second] = false;
0246           }
0247         }
0248       }
0249     }
0250   }
0251 }