File indexing completed on 2024-09-07 04:37:50
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
0039
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
0061 constexpr double delta_cPhi = 0.00638;
0062
0063 double defi = 0;
0064 if (eta < 0)
0065 phi += delta_cPhi;
0066 if (phi >= -pi && phi <= pi) {
0067
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 {
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 }
0083 }
0084 } else {
0085 defi = 0;
0086 }
0087
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 }
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
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
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;
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
0156
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
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
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
0198
0199 const double rhenergy = rechit.energy();
0200
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();
0208 surroundingEnergy += sum;
0209 }
0210
0211
0212
0213 const double fraction1 = surroundingEnergy / rhenergy;
0214
0215
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 }
0228
0229 if (mask[i] && rhenergy > clean._doubleSpikeThresh) {
0230
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
0247 if (enmax > 0.0) {
0248 double surroundingEnergyj = 0.0;
0249 auto const& neighbours4j = hits[mostEnergeticNeighbour].neighbours4();
0250 for (auto k : neighbours4j) {
0251
0252 surroundingEnergyj += hits[k].energy();
0253 }
0254
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 }
0271 }
0272 }
0273 }