JetMET/src/HLTHcalTowerNoiseCleaner.cc

0001 // -*- C++ -*-
0002 //
0003 // Class:      HLTHcalTowerNoiseCleaner
0004 //
0005 /**\class HLTHcalTowerNoiseCleaner
0006
0007  Description: HLT filter module for cleaning HCal Noise from MET or MHT
0008
0009  Implementation:
0010      <Notes on implementation>
0011 */
0012 //
0013 // Original Author:  Alexander Mott
0014 //         Created:  Mon Nov 21 11:32:00 CEST 2011
0015 //
0016 //
0017 //
0018 #include <set>
0019
0020 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0021 #include "FWCore/Framework/interface/Event.h"
0022 #include "FWCore/Framework/interface/EventSetup.h"
0023 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0024 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0025 #include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
0026 #include "FWCore/Utilities/interface/InputTag.h"
0027 #include "DataFormats/Common/interface/Handle.h"
0028 #include "DataFormats/Candidate/interface/Candidate.h"
0029 #include "DataFormats/JetReco/interface/Jet.h"
0030 #include "DataFormats/JetReco/interface/CaloJet.h"
0031 #include "DataFormats/CaloTowers/interface/CaloTowerDetId.h"
0032 #include "DataFormats/CaloTowers/interface/CaloTowerCollection.h"
0033 #include "HLTrigger/JetMET/interface/HLTHcalTowerNoiseCleaner.h"
0034
0035 HLTHcalTowerNoiseCleaner::HLTHcalTowerNoiseCleaner(const edm::ParameterSet& iConfig)
0036     : hcalRecNumberingRecordToken_(esConsumes()),
0037       HcalNoiseRBXCollectionTag_(iConfig.getParameter<edm::InputTag>("HcalNoiseRBXCollection")),
0038       TowerCollectionTag_(iConfig.getParameter<edm::InputTag>("CaloTowerCollection")),
0039       severity_(iConfig.getParameter<int>("severity")),
0040       maxNumRBXs_(iConfig.getParameter<int>("maxNumRBXs")),
0041       numRBXsToConsider_(iConfig.getParameter<int>("numRBXsToConsider")),
0042       needEMFCoincidence_(iConfig.getParameter<bool>("needEMFCoincidence")),
0043       minRBXEnergy_(iConfig.getParameter<double>("minRBXEnergy")),
0044       minRatio_(iConfig.getParameter<double>("minRatio")),
0045       maxRatio_(iConfig.getParameter<double>("maxRatio")),
0046       minHPDHits_(iConfig.getParameter<int>("minHPDHits")),
0047       minRBXHits_(iConfig.getParameter<int>("minRBXHits")),
0048       minHPDNoOtherHits_(iConfig.getParameter<int>("minHPDNoOtherHits")),
0049       minZeros_(iConfig.getParameter<int>("minZeros")),
0050       minHighEHitTime_(iConfig.getParameter<double>("minHighEHitTime")),
0051       maxHighEHitTime_(iConfig.getParameter<double>("maxHighEHitTime")),
0052       maxRBXEMF_(iConfig.getParameter<double>("maxRBXEMF")),
0053       minRecHitE_(iConfig.getParameter<double>("minRecHitE")),
0054       minLowHitE_(iConfig.getParameter<double>("minLowHitE")),
0055       minHighHitE_(iConfig.getParameter<double>("minHighHitE")),
0056       minR45HitE_(5.0),
0057       TS4TS5EnergyThreshold_(iConfig.getParameter<double>("TS4TS5EnergyThreshold")) {
0058   std::vector<double> TS4TS5UpperThresholdTemp = iConfig.getParameter<std::vector<double> >("TS4TS5UpperThreshold");
0059   std::vector<double> TS4TS5UpperCutTemp = iConfig.getParameter<std::vector<double> >("TS4TS5UpperCut");
0060   std::vector<double> TS4TS5LowerThresholdTemp = iConfig.getParameter<std::vector<double> >("TS4TS5LowerThreshold");
0061   std::vector<double> TS4TS5LowerCutTemp = iConfig.getParameter<std::vector<double> >("TS4TS5LowerCut");
0062
0063   for (int i = 0; i < (int)TS4TS5UpperThresholdTemp.size() && i < (int)TS4TS5UpperCutTemp.size(); i++)
0064     TS4TS5UpperCut_.push_back(std::pair<double, double>(TS4TS5UpperThresholdTemp[i], TS4TS5UpperCutTemp[i]));
0065   sort(TS4TS5UpperCut_.begin(), TS4TS5UpperCut_.end());
0066
0067   for (int i = 0; i < (int)TS4TS5LowerThresholdTemp.size() && i < (int)TS4TS5LowerCutTemp.size(); i++)
0068     TS4TS5LowerCut_.push_back(std::pair<double, double>(TS4TS5LowerThresholdTemp[i], TS4TS5LowerCutTemp[i]));
0069   sort(TS4TS5LowerCut_.begin(), TS4TS5LowerCut_.end());
0070
0071   m_theHcalNoiseToken = consumes<reco::HcalNoiseRBXCollection>(HcalNoiseRBXCollectionTag_);
0072   m_theCaloTowerCollectionToken = consumes<CaloTowerCollection>(TowerCollectionTag_);
0073
0074   if (iConfig.existsAs<double>("minR45HitE"))
0075     minR45HitE_ = iConfig.getParameter<double>("minR45HitE");
0076
0077   produces<CaloTowerCollection>();
0078 }
0079
0080 HLTHcalTowerNoiseCleaner::~HLTHcalTowerNoiseCleaner() = default;
0081
0082 void HLTHcalTowerNoiseCleaner::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0083   edm::ParameterSetDescription desc;
0084   desc.add<edm::InputTag>("HcalNoiseRBXCollection", edm::InputTag("hltHcalNoiseInfoProducer"));
0085   desc.add<edm::InputTag>("CaloTowerCollection", edm::InputTag("hltTowerMakerForAll"));
0086   desc.add<double>("maxTowerNoiseEnergyFraction", 0.5);
0087   desc.add<int>("severity", 1);
0088   desc.add<int>("maxNumRBXs", 2);
0089   desc.add<int>("numRBXsToConsider", 2);
0090   desc.add<bool>("needEMFCoincidence", true);
0091   desc.add<double>("minRBXEnergy", 50.0);
0092   desc.add<double>("minRatio", -999.);
0093   desc.add<double>("maxRatio", 999.);
0094   desc.add<int>("minHPDHits", 17);
0095   desc.add<int>("minRBXHits", 999);
0096   desc.add<int>("minHPDNoOtherHits", 10);
0097   desc.add<int>("minZeros", 10);
0098   desc.add<double>("minHighEHitTime", -9999.0);
0099   desc.add<double>("maxHighEHitTime", 9999.0);
0100   desc.add<double>("maxRBXEMF", 0.02);
0101   desc.add<double>("minRecHitE", 1.5);
0102   desc.add<double>("minLowHitE", 10.0);
0103   desc.add<double>("minHighHitE", 25.0);
0104   desc.add<double>("minR45HitE", 5.0);
0105   desc.add<double>("TS4TS5EnergyThreshold", 50.0);
0106
0107   double TS4TS5UpperThresholdArray[5] = {70, 90, 100, 400, 4000};
0108   double TS4TS5UpperCutArray[5] = {1, 0.8, 0.75, 0.72, 0.72};
0109   double TS4TS5LowerThresholdArray[7] = {100, 120, 150, 200, 300, 400, 500};
0110   double TS4TS5LowerCutArray[7] = {-1, -0.7, -0.4, -0.2, -0.08, 0, 0.1};
0111   std::vector<double> TS4TS5UpperThreshold(TS4TS5UpperThresholdArray, TS4TS5UpperThresholdArray + 5);
0112   std::vector<double> TS4TS5UpperCut(TS4TS5UpperCutArray, TS4TS5UpperCutArray + 5);
0113   std::vector<double> TS4TS5LowerThreshold(TS4TS5LowerThresholdArray, TS4TS5LowerThresholdArray + 7);
0114   std::vector<double> TS4TS5LowerCut(TS4TS5LowerCutArray, TS4TS5LowerCutArray + 7);
0115
0116   desc.add<std::vector<double> >("TS4TS5UpperThreshold", TS4TS5UpperThreshold);
0117   desc.add<std::vector<double> >("TS4TS5UpperCut", TS4TS5UpperCut);
0118   desc.add<std::vector<double> >("TS4TS5LowerThreshold", TS4TS5LowerThreshold);
0119   desc.add<std::vector<double> >("TS4TS5LowerCut", TS4TS5LowerCut);
0120   descriptions.add("hltHcalTowerNoiseCleaner", desc);
0121 }
0122
0123 //
0124 // member functions
0125 //
0126
0127 void HLTHcalTowerNoiseCleaner::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0128   using namespace reco;
0129
0130   //get the calo MET / MHT
0131   edm::Handle<CaloTowerCollection> tower_h;
0132   iEvent.getByToken(m_theCaloTowerCollectionToken, tower_h);
0133
0134   std::set<unsigned int> noisyTowers;
0135
0136   if (not tower_h.isValid()) {  //No towers MET, don't do anything and accept the event
0137     edm::LogError("HLTHcalTowerNoiseCleaner") << "Input Tower Collection is not Valid";
0138     return;
0139   }
0140
0141   //get the calotower topology
0142   auto const& caloTowerTopology = iSetup.getHandle(hcalRecNumberingRecordToken_);
0143
0144   // get the RBXs produced by RecoMET/METProducers/HcalNoiseInfoProducer
0145   edm::Handle<HcalNoiseRBXCollection> rbxs_h;
0146   iEvent.getByToken(m_theHcalNoiseToken, rbxs_h);
0147   if (!rbxs_h.isValid()) {
0148     edm::LogWarning("HLTHcalTowerNoiseCleaner")
0149         << "Could not find HcalNoiseRBXCollection product named " << HcalNoiseRBXCollectionTag_ << "." << std::endl;
0150     severity_ = 0;
0151   }
0152
0153   // create a sorted set of the RBXs, ordered by energy
0154   noisedataset_t data;
0155   for (auto const& rbx : *rbxs_h) {
0156     CommonHcalNoiseRBXData d(rbx,
0157                              minRecHitE_,
0158                              minLowHitE_,
0159                              minHighHitE_,
0160                              TS4TS5EnergyThreshold_,
0161                              TS4TS5UpperCut_,
0162                              TS4TS5LowerCut_,
0163                              minR45HitE_);
0164     data.insert(d);
0165   }
0166
0167   // data is now sorted by RBX energy
0168   // only consider top N=numRBXsToConsider_ energy RBXs
0169   if (severity_ > 0) {
0170     for (auto const& it : data) {
0171       bool passFilter = true;
0172       bool passEMF = true;
0173       if (it.energy() > minRBXEnergy_) {
0174         if (it.validRatio() && it.ratio() < minRatio_)
0175           passFilter = false;
0176         else if (it.validRatio() && it.ratio() > maxRatio_)
0177           passFilter = false;
0178         else if (it.numHPDHits() >= minHPDHits_)
0179           passFilter = false;
0180         else if (it.numRBXHits() >= minRBXHits_)
0181           passFilter = false;
0182         else if (it.numHPDNoOtherHits() >= minHPDNoOtherHits_)
0183           passFilter = false;
0184         else if (it.numZeros() >= minZeros_)
0185           passFilter = false;
0186         else if (it.minHighEHitTime() < minHighEHitTime_)
0187           passFilter = false;
0188         else if (it.maxHighEHitTime() > maxHighEHitTime_)
0189           passFilter = false;
0190         else if (!it.PassTS4TS5())
0191           passFilter = false;
0192
0193         if (it.RBXEMF() < maxRBXEMF_) {
0194           passEMF = false;
0195         }
0196       }
0197
0198       if ((needEMFCoincidence_ && !passEMF && !passFilter) ||
0199           (!needEMFCoincidence_ && !passFilter)) {  // check for noise
0200         LogDebug("") << "HLTHcalTowerNoiseCleaner debug: Found a noisy RBX: "
0201                      << "energy=" << it.energy() << "; "
0202                      << "ratio=" << it.ratio() << "; "
0203                      << "# RBX hits=" << it.numRBXHits() << "; "
0204                      << "# HPD hits=" << it.numHPDHits() << "; "
0205                      << "# Zeros=" << it.numZeros() << "; "
0206                      << "min time=" << it.minHighEHitTime() << "; "
0207                      << "max time=" << it.maxHighEHitTime() << "; "
0208                      << "passTS4TS5=" << it.PassTS4TS5() << "; "
0209                      << "RBX EMF=" << it.RBXEMF() << std::endl;
0210         // add calotowers associated with this RBX to the noise list
0211         edm::RefVector<CaloTowerCollection> noiseTowers = it.rbxTowers();
0212         edm::RefVector<CaloTowerCollection>::const_iterator noiseTowersIt;
0213         //add these calotowers to the noisy list
0214         for (noiseTowersIt = noiseTowers.begin(); noiseTowersIt != noiseTowers.end(); noiseTowersIt++) {
0215           edm::Ref<edm::SortedCollection<CaloTower> > tower_ref = *noiseTowersIt;
0216           CaloTowerDetId id = tower_ref->id();
0217           noisyTowers.insert(caloTowerTopology->denseIndex(id));
0218         }
0219       }
0220     }  // done with noise loop
0221   }    //if(severity_>0)
0222
0223   //output collection
0224   std::unique_ptr<CaloTowerCollection> OutputTowers(new CaloTowerCollection());
0225
0226   CaloTowerCollection::const_iterator inTowersIt;
0227
0228   for (inTowersIt = tower_h->begin(); inTowersIt != tower_h->end(); inTowersIt++) {
0229     const CaloTower& tower = (*inTowersIt);
0230     CaloTowerDetId id = tower.id();
0231     if (noisyTowers.find(caloTowerTopology->denseIndex(id)) == noisyTowers.end()) {  // the tower is not noisy
0232       OutputTowers->push_back(*inTowersIt);
0233     }
0234   }
0235   iEvent.put(std::move(OutputTowers));
0236 }