PFClusterProducer/plugins/SpikeAndDoubleSpikeCleaner.cc

0001 #include "DataFormats/HcalRecHit/interface/HFRecHit.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 SpikeAndDoubleSpikeCleaner : public RecHitTopologicalCleanerBase {
0009 public:
0010   struct spike_cleaning {
0011     double _singleSpikeThresh;
0012     double _minS4S1_a;
0013     double _minS4S1_b;
0014     double _doubleSpikeS6S2;
0015     double _eneThreshMod;
0016     double _fracThreshMod;
0017     double _doubleSpikeThresh;
0018   };
0019
0020   SpikeAndDoubleSpikeCleaner(const edm::ParameterSet& conf, edm::ConsumesCollector& cc);
0021   SpikeAndDoubleSpikeCleaner(const SpikeAndDoubleSpikeCleaner&) = delete;
0022   SpikeAndDoubleSpikeCleaner& operator=(const SpikeAndDoubleSpikeCleaner&) = delete;
0023
0024   void clean(const edm::Handle<reco::PFRecHitCollection>& input, std::vector<bool>& mask) override;
0025
0026 private:
0027   const std::unordered_map<std::string, int> _layerMap;
0028   std::unordered_map<int, spike_cleaning> _thresholds;
0029 };
0030
0031 DEFINE_EDM_PLUGIN(RecHitTopologicalCleanerFactory, SpikeAndDoubleSpikeCleaner, "SpikeAndDoubleSpikeCleaner");
0032
0033 namespace {
0034   std::pair<double, double> dCrack(double phi, double eta) {
0035     constexpr double oneOverCrystalSize = 1.0 / 0.0175;
0036     constexpr double pi = M_PI;
0037     constexpr double twopi = 2 * pi;
0038     // the result below is from unrolling
0039     // the lazy-eval loop in PFClusterAlgo::dCrack
0040     constexpr double cPhi[18] = {2.97025,
0041                                  2.621184149601134,
0042                                  2.272118299202268,
0043                                  1.9230524488034024,
0044                                  1.5739865984045365,
0045                                  1.2249207480056705,
0046                                  0.8758548976068048,
0047                                  0.5267890472079388,
0048                                  0.1777231968090729,
0049                                  -0.17134265358979306,
0050                                  -0.520408503988659,
0051                                  -0.8694743543875245,
0052                                  -1.2185402047863905,
0053                                  -1.5676060551852569,
0054                                  -1.9166719055841224,
0055                                  -2.265737755982988,
0056                                  -2.6148036063818543,
0057                                  -2.9638694567807207};
0058     constexpr double cEta[9] = {
0059         0.0, 4.44747e-01, -4.44747e-01, 7.92824e-01, -7.92824e-01, 1.14090e+00, -1.14090e+00, 1.47464e+00, -1.47464e+00};
0060     // shift for eta < 0
0061     constexpr double delta_cPhi = 0.00638;
0062     // let's calculate dphi
0063     double defi = 0;
0064     if (eta < 0)
0065       phi += delta_cPhi;
0066     if (phi >= -pi && phi <= pi) {
0067       //the problem of the extrema
0068       if (phi < cPhi[17] || phi >= cPhi[0]) {
0069         if (phi < 0)
0070           phi += 2 * pi;
0071         defi = std::min(std::abs(phi - cPhi[0]), std::abs(phi - cPhi[17] - twopi));
0072       } else {  // between these extrema
0073         bool OK = false;
0074         unsigned i = 16;
0075         while (!OK) {
0076           if (phi < cPhi[i]) {
0077             defi = std::min(std::abs(phi - cPhi[i + 1]), std::abs(phi - cPhi[i]));
0078             OK = true;
0079           } else {
0080             i -= 1;
0081           }
0082         }  // end while
0083       }
0084     } else {  // if there's a problem assume we're in a crack
0085       defi = 0;
0086     }
0087     // let's calculate deta
0088     double deta = 999.0;
0089     for (const double etaGap : cEta) {
0090       deta = std::min(deta, std::abs(eta - etaGap));
0091     }
0092     defi *= oneOverCrystalSize;
0093     deta *= oneOverCrystalSize;
0094     return std::make_pair(defi, deta);
0095   }
0096 }  // namespace
0097
0098 SpikeAndDoubleSpikeCleaner::SpikeAndDoubleSpikeCleaner(const edm::ParameterSet& conf, edm::ConsumesCollector& cc)
0099     : RecHitTopologicalCleanerBase(conf, cc),
0100       _layerMap({{"PS2", (int)PFLayer::PS2},
0101                  {"PS1", (int)PFLayer::PS1},
0102                  {"ECAL_ENDCAP", (int)PFLayer::ECAL_ENDCAP},
0103                  {"ECAL_BARREL", (int)PFLayer::ECAL_BARREL},
0104                  {"NONE", (int)PFLayer::NONE},
0105                  {"HCAL_BARREL1", (int)PFLayer::HCAL_BARREL1},
0106                  {"HCAL_BARREL2_RING0", (int)PFLayer::HCAL_BARREL2},
0107                  // hack to deal with ring1 in HO
0108                  {"HCAL_BARREL2_RING1", 100 * (int)PFLayer::HCAL_BARREL2},
0109                  {"HCAL_ENDCAP", (int)PFLayer::HCAL_ENDCAP},
0110                  {"HF_EM", (int)PFLayer::HF_EM},
0111                  {"HF_HAD", (int)PFLayer::HF_HAD}}) {
0112   const std::vector<edm::ParameterSet>& thresholds = conf.getParameterSetVector("cleaningByDetector");
0113   for (const auto& pset : thresholds) {
0114     spike_cleaning info;
0115     const std::string& det = pset.getParameter<std::string>("detector");
0116     info._minS4S1_a = pset.getParameter<double>("minS4S1_a");
0117     info._minS4S1_b = pset.getParameter<double>("minS4S1_b");
0118     info._doubleSpikeS6S2 = pset.getParameter<double>("doubleSpikeS6S2");
0119     info._eneThreshMod = pset.getParameter<double>("energyThresholdModifier");
0120     info._fracThreshMod = pset.getParameter<double>("fractionThresholdModifier");
0121     info._doubleSpikeThresh = pset.getParameter<double>("doubleSpikeThresh");
0122     info._singleSpikeThresh = pset.getParameter<double>("singleSpikeThresh");
0123     auto entry = _layerMap.find(det);
0124     if (entry == _layerMap.end()) {
0125       throw cms::Exception("InvalidDetectorLayer") << "Detector layer : " << det << " is not in the list of recognized"
0126                                                    << " detector layers!";
0127     }
0128     _thresholds.emplace(_layerMap.find(det)->second, info);
0129   }
0130 }
0131
0132 void SpikeAndDoubleSpikeCleaner::clean(const edm::Handle<reco::PFRecHitCollection>& input, std::vector<bool>& mask) {
0133   //need to run over energy sorted rechits
0134   auto const& hits = *input;
0135   std::vector<unsigned> ordered_hits(hits.size());
0136   for (unsigned i = 0; i < hits.size(); ++i)
0137     ordered_hits[i] = i;
0138   std::sort(ordered_hits.begin(), ordered_hits.end(), [&](unsigned i, unsigned j) {
0139     return hits[i].energy() > hits[j].energy();
0140   });
0141
0142   for (const auto& idx : ordered_hits) {
0143     const unsigned i = idx;
0144     if (!mask[i])
0145       continue;  // don't need to re-mask things :-)
0146     const reco::PFRecHit& rechit = hits[i];
0147     int hitlayer = (int)rechit.layer();
0148     if (hitlayer == PFLayer::HCAL_BARREL2 && std::abs(rechit.positionREP().eta()) > 0.34) {
0149       hitlayer *= 100;
0150     }
0151     const spike_cleaning& clean = _thresholds.find(hitlayer)->second;
0152     if (rechit.energy() < clean._singleSpikeThresh)
0153       continue;
0154
0155     //Fix needed for HF.  Here, we find the (up to) five companion rechits
0156     //to work in conjunction with the neighbours4() implementation below for the full HF surrounding energy
0157     float compsumE = 0.0;
0158     if ((hitlayer == PFLayer::HF_EM || hitlayer == PFLayer::HF_HAD)) {
0159       int comp = 1;
0160       if (hitlayer == PFLayer::HF_EM)
0161         comp = 2;
0162       const HcalDetId& detid = (HcalDetId)rechit.detId();
0163       int heta = detid.ieta();
0164       int hphi = detid.iphi();
0165
0166       //At eta>39, phi segmentation changes
0167       int predphi = 2;
0168       if (std::abs(heta) > 39)
0169         predphi = 4;
0170
0171       int curphiL = hphi - predphi;
0172       int curphiH = hphi + predphi;
0173
0174       //HcalDetId valid phi range (0-72)
0175       while (curphiL < 0)
0176         curphiL += 72;
0177       while (curphiH > 72)
0178         curphiH -= 72;
0179
0180       std::pair<std::vector<int>, std::vector<int>> phietas({heta, heta + 1, heta - 1, heta, heta},
0181                                                             {hphi, hphi, hphi, curphiL, curphiH});
0182
0183       std::vector<uint32_t> rawDetIds;
0184       for (unsigned in = 0; in < phietas.first.size(); in++) {
0185         HcalDetId tempID(HcalForward, phietas.first[in], phietas.second[in], comp);
0186         rawDetIds.push_back(tempID.rawId());
0187       }
0188
0189       for (const auto& jdx : ordered_hits) {
0190         const unsigned j = jdx;
0191         const reco::PFRecHit& matchrechit = hits[j];
0192         for (const auto& iID : rawDetIds)
0193           if (iID == matchrechit.detId())
0194             compsumE += matchrechit.energy();
0195       }
0196     }
0197     //End of fix needed for HF
0198
0199     const double rhenergy = rechit.energy();
0200     // single spike cleaning
0201     auto const& neighbours4 = rechit.neighbours4();
0202     double surroundingEnergy = compsumE;
0203     for (auto k : neighbours4) {
0204       if (!mask[k])
0205         continue;
0206       auto const& neighbour = hits[k];
0207       const double sum = neighbour.energy();  //energyUp is just rechit energy?
0208       surroundingEnergy += sum;
0209     }
0210     //   wannaBeSeed.energyUp()/wannaBeSeed.energy() : 1.;
0211     // Fraction 1 is the balance between the hit and its neighbours
0212     // from both layers
0213     const double fraction1 = surroundingEnergy / rhenergy;
0214     // removed spurious comments from old pfcluster algo...
0215     // look there if you want more history
0216     const double f1Cut = (clean._minS4S1_a * std::log10(rechit.energy()) + clean._minS4S1_b);
0217     if (fraction1 < f1Cut) {
0218       const double eta = rechit.positionREP().eta();
0219       const double aeta = std::abs(eta);
0220       const double phi = rechit.positionREP().phi();
0221       std::pair<double, double> dcr = dCrack(phi, eta);
0222       const double dcrmin = (rechit.layer() == PFLayer::ECAL_BARREL ? std::min(dcr.first, dcr.second) : dcr.second);
0223       if (aeta < 5.0 && ((aeta < 2.85 && dcrmin > 1.0) || (rhenergy > clean._eneThreshMod * clean._singleSpikeThresh &&
0224                                                            fraction1 < f1Cut / clean._fracThreshMod))) {
0225         mask[i] = false;
0226       }
0227     }  //if initial fraction cut (single spike)
0228     // double spike removal
0229     if (mask[i] && rhenergy > clean._doubleSpikeThresh) {
0230       //Determine energy surrounding the seed and the most energetic neighbour
0231       double surroundingEnergyi = 0.0;
0232       double enmax = -999.0;
0233       unsigned int mostEnergeticNeighbour = 0;
0234       auto const& neighbours4i = rechit.neighbours4();
0235       for (auto k : neighbours4i) {
0236         if (!mask[k])
0237           continue;
0238         auto const& neighbour = hits[k];
0239         const double nenergy = neighbour.energy();
0240         surroundingEnergyi += nenergy;
0241         if (nenergy > enmax) {
0242           enmax = nenergy;
0243           mostEnergeticNeighbour = k;
0244         }
0245       }
0246       // is there an energetic neighbour
0247       if (enmax > 0.0) {
0248         double surroundingEnergyj = 0.0;
0249         auto const& neighbours4j = hits[mostEnergeticNeighbour].neighbours4();
0250         for (auto k : neighbours4j) {
0251           //if( !mask[k] &&  k != i) continue; // leave out?
0252           surroundingEnergyj += hits[k].energy();
0253         }
0254         // The energy surrounding the double spike candidate
0255         const double surroundingEnergyFraction =
0256             (surroundingEnergyi + surroundingEnergyj) / (rechit.energy() + hits[mostEnergeticNeighbour].energy()) - 1.;
0257         if (surroundingEnergyFraction < clean._doubleSpikeS6S2) {
0258           const double eta = rechit.positionREP().eta();
0259           const double aeta = std::abs(eta);
0260           const double phi = rechit.positionREP().phi();
0261           std::pair<double, double> dcr = dCrack(phi, eta);
0262           const double dcrmin = (rechit.layer() == PFLayer::ECAL_BARREL ? std::min(dcr.first, dcr.second) : dcr.second);
0263           if (aeta < 5.0 && ((aeta < 2.85 && dcrmin > 1.0) ||
0264                              (rhenergy > clean._eneThreshMod * clean._doubleSpikeThresh &&
0265                               surroundingEnergyFraction < clean._doubleSpikeS6S2 / clean._fracThreshMod))) {
0266             mask[i] = false;
0267             mask[mostEnergeticNeighbour] = false;
0268           }
0269         }
0270       }  // was there an energetic neighbour ?
0271     }    // if double spike thresh
0272   }      // rechit loop
0273 }