Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-03-25 23:59:52

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