Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 14:20:31

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