File indexing completed on 2024-04-06 12:18:33
0001
0002
0003
0004
0005
0006
0007 #include <array>
0008 #include <cmath>
0009 #include <utility>
0010
0011 #include "DataFormats/Common/interface/Handle.h"
0012 #include "DataFormats/HcalDetId/interface/HcalDetId.h"
0013 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0014 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0015 #include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
0016 #include "FWCore/Utilities/interface/InputTag.h"
0017 #include "HLTrigger/JetMET/interface/HLTHPDFilter.h"
0018
0019 namespace {
0020 enum Partition { HBM = 0, HBP = 1, HEM = 2, HEP = 3 };
0021 std::pair<Partition, int> hpdId(HcalDetId fId) {
0022 int hpd = fId.iphi();
0023 Partition partition = HBM;
0024 if (fId.subdet() == HcalBarrel) {
0025 partition = fId.ieta() > 0 ? HBP : HBM;
0026 } else if (fId.subdet() == HcalEndcap) {
0027 partition = fId.ieta() > 0 ? HEP : HEM;
0028 if ((fId.iphi() - 1) % 4 < 2) {
0029 switch (fId.ieta()) {
0030 case 22:
0031 case 24:
0032 case 26:
0033 case 28:
0034 hpd = +1;
0035 break;
0036 case 29:
0037 if (fId.depth() == 1 || fId.depth() == 3)
0038 hpd += 1;
0039 break;
0040 default:
0041 break;
0042 }
0043 } else {
0044 switch (fId.ieta()) {
0045 case 21:
0046 case 23:
0047 case 25:
0048 case 27:
0049 hpd += 1;
0050 break;
0051 case 29:
0052 if (fId.depth() == 2)
0053 hpd += 1;
0054 break;
0055 default:
0056 break;
0057 }
0058 }
0059 }
0060 return std::pair<Partition, int>(partition, hpd);
0061 }
0062 }
0063
0064 HLTHPDFilter::HLTHPDFilter(const edm::ParameterSet& iConfig)
0065 : mInputTag(iConfig.getParameter<edm::InputTag>("inputTag")),
0066 mEnergyThreshold(iConfig.getParameter<double>("energy")),
0067 mHPDSpikeEnergyThreshold(iConfig.getParameter<double>("hpdSpikeEnergy")),
0068 mHPDSpikeIsolationEnergyThreshold(iConfig.getParameter<double>("hpdSpikeIsolationEnergy")),
0069 mRBXSpikeEnergyThreshold(iConfig.getParameter<double>("rbxSpikeEnergy")),
0070 mRBXSpikeUnbalanceThreshold(iConfig.getParameter<double>("rbxSpikeUnbalance")) {
0071 m_theRecHitCollectionToken = consumes<HBHERecHitCollection>(mInputTag);
0072 }
0073
0074 HLTHPDFilter::~HLTHPDFilter() = default;
0075
0076 void HLTHPDFilter::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0077 edm::ParameterSetDescription desc;
0078 desc.add<edm::InputTag>("inputTag", edm::InputTag("hltHbhereco"));
0079 desc.add<double>("energy", -99.0);
0080 desc.add<double>("hpdSpikeEnergy", 10.0);
0081 desc.add<double>("hpdSpikeIsolationEnergy", 1.0);
0082 desc.add<double>("rbxSpikeEnergy", 50.0);
0083 desc.add<double>("rbxSpikeUnbalance", 0.2);
0084 descriptions.add("hltHPDFilter", desc);
0085 }
0086
0087 bool HLTHPDFilter::filter(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0088 if (mHPDSpikeEnergyThreshold <= 0 && mRBXSpikeEnergyThreshold <= 0)
0089 return true;
0090
0091 edm::Handle<HBHERecHitCollection> hbhe;
0092 iEvent.getByToken(m_theRecHitCollectionToken, hbhe);
0093
0094
0095 std::array<std::array<float, 73>, 4> hpdEnergy;
0096 for (auto& i : hpdEnergy)
0097 for (size_t j = 0; j < 73; ++j)
0098 i[j] = 0;
0099
0100
0101 for (unsigned i = 0; i < hbhe->size(); ++i) {
0102 if ((*hbhe)[i].energy() > mEnergyThreshold) {
0103 std::pair<Partition, int> hpd = hpdId((*hbhe)[i].id());
0104 hpdEnergy[int(hpd.first)][hpd.second] += (*hbhe)[i].energy();
0105 }
0106 }
0107
0108
0109 if (mHPDSpikeEnergyThreshold > 0) {
0110 for (auto& partition : hpdEnergy) {
0111 for (size_t i = 1; i < 73; ++i) {
0112 if (partition[i] > mHPDSpikeEnergyThreshold) {
0113 int hpdPlus = i + 1;
0114 if (hpdPlus == 73)
0115 hpdPlus = 1;
0116 int hpdMinus = i - 1;
0117 if (hpdMinus == 0)
0118 hpdMinus = 72;
0119 double maxNeighborEnergy = fmax(partition[hpdPlus], partition[hpdMinus]);
0120 if (maxNeighborEnergy < mHPDSpikeIsolationEnergyThreshold)
0121 return false;
0122 }
0123 }
0124 }
0125 }
0126
0127
0128 if (mRBXSpikeEnergyThreshold > 0) {
0129 for (auto& partition : hpdEnergy) {
0130 for (size_t rbx = 1; rbx < 19; ++rbx) {
0131 int ifirst = (rbx - 1) * 4 - 1;
0132 int iend = (rbx - 1) * 4 + 3;
0133 double minEnergy = 0;
0134 double maxEnergy = -1;
0135 double totalEnergy = 0;
0136 for (int irm = ifirst; irm < iend; ++irm) {
0137 int hpd = irm;
0138 if (hpd <= 0)
0139 hpd = 72 + hpd;
0140 totalEnergy += partition[hpd];
0141 if (minEnergy > maxEnergy) {
0142 minEnergy = maxEnergy = partition[hpd];
0143 } else {
0144 if (partition[hpd] < minEnergy)
0145 minEnergy = partition[hpd];
0146 if (partition[hpd] > maxEnergy)
0147 maxEnergy = partition[hpd];
0148 }
0149 }
0150 if (totalEnergy > mRBXSpikeEnergyThreshold) {
0151 if (minEnergy / maxEnergy > mRBXSpikeUnbalanceThreshold)
0152 return false;
0153 }
0154 }
0155 }
0156 }
0157 return true;
0158 }