JetMET/src/HLTHPDFilter.cc

0001 /* \class HLTHPDFilter
0002  *
0003  *
0004  * Fedor Ratnikov (UMd) May 19, 2008
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) {  // 1,2
0029         switch (fId.ieta()) {          // 1->2
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 {                 // 3,4
0044         switch (fId.ieta()) {  // 3->4
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 }  // namespace
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;  // nothing to filter
0090   // get hits
0091   edm::Handle<HBHERecHitCollection> hbhe;
0092   iEvent.getByToken(m_theRecHitCollectionToken, hbhe);
0093
0094   // collect energies
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   // select hist above threshold
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   // not single HPD spike
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;  // HPD spike found
0122         }
0123       }
0124     }
0125   }
0126
0127   // not RBX flash
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;  // likely HPF flash
0153         }
0154       }
0155     }
0156   }
0157   return true;
0158 }