Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2023-03-17 11:09:34 

0001 // -*- C++ -*-
0002 //
0003 // Class:      HLTHcalMETNoiseCleaner
0004 //
0005 /**\class HLTHcalMETNoiseCleaner
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 
0019 #include "HLTrigger/JetMET/interface/HLTHcalMETNoiseCleaner.h"
0020 
0021 #include "DataFormats/Common/interface/Handle.h"
0022 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0023 #include "FWCore/Framework/interface/ESHandle.h"
0024 #include "FWCore/Framework/interface/Event.h"
0025 #include "FWCore/Framework/interface/EventSetup.h"
0026 
0027 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0028 
0029 #include "DataFormats/Candidate/interface/Candidate.h"
0030 #include "DataFormats/METReco/interface/SpecificCaloMETData.h"
0031 #include "DataFormats/CaloTowers/interface/CaloTower.h"
0032 #include "DataFormats/CaloTowers/interface/CaloTowerCollection.h"
0033 #include "DataFormats/Math/interface/LorentzVector.h"
0034 #include "DataFormats/Math/interface/Point3D.h"
0035 
0036 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0037 #include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
0038 #include "FWCore/Utilities/interface/InputTag.h"
0039 
0040 #include <iostream>
0041 #include <string>
0042 #include <fstream>
0043 #include <TVector3.h>
0044 #include <TLorentzVector.h>
0045 //#include <Point.h>
0046 
0047 HLTHcalMETNoiseCleaner::HLTHcalMETNoiseCleaner(const edm::ParameterSet& iConfig)
0048     : HcalNoiseRBXCollectionTag_(iConfig.getParameter<edm::InputTag>("HcalNoiseRBXCollection")),
0049       CaloMetCollectionTag_(iConfig.getParameter<edm::InputTag>("CaloMetCollection")),
0050       CaloMetCut_(iConfig.getParameter<double>("CaloMetCut")),
0051       severity_(iConfig.getParameter<int>("severity")),
0052       maxNumRBXs_(iConfig.getParameter<int>("maxNumRBXs")),
0053       numRBXsToConsider_(iConfig.getParameter<int>("numRBXsToConsider")),
0054       accept2NoiseRBXEvents_(iConfig.getParameter<bool>("accept2NoiseRBXEvents")),
0055       needEMFCoincidence_(iConfig.getParameter<bool>("needEMFCoincidence")),
0056       minRBXEnergy_(iConfig.getParameter<double>("minRBXEnergy")),
0057       minRatio_(iConfig.getParameter<double>("minRatio")),
0058       maxRatio_(iConfig.getParameter<double>("maxRatio")),
0059       minHPDHits_(iConfig.getParameter<int>("minHPDHits")),
0060       minRBXHits_(iConfig.getParameter<int>("minRBXHits")),
0061       minHPDNoOtherHits_(iConfig.getParameter<int>("minHPDNoOtherHits")),
0062       minZeros_(iConfig.getParameter<int>("minZeros")),
0063       minHighEHitTime_(iConfig.getParameter<double>("minHighEHitTime")),
0064       maxHighEHitTime_(iConfig.getParameter<double>("maxHighEHitTime")),
0065       maxRBXEMF_(iConfig.getParameter<double>("maxRBXEMF")),
0066       minRecHitE_(iConfig.getParameter<double>("minRecHitE")),
0067       minLowHitE_(iConfig.getParameter<double>("minLowHitE")),
0068       minHighHitE_(iConfig.getParameter<double>("minHighHitE")),
0069       minR45HitE_(iConfig.getParameter<double>("minR45HitE")),
0070       TS4TS5EnergyThreshold_(iConfig.getParameter<double>("TS4TS5EnergyThreshold")) {
0071   std::vector<double> TS4TS5UpperThresholdTemp = iConfig.getParameter<std::vector<double> >("TS4TS5UpperThreshold");
0072   std::vector<double> TS4TS5UpperCutTemp = iConfig.getParameter<std::vector<double> >("TS4TS5UpperCut");
0073   std::vector<double> TS4TS5LowerThresholdTemp = iConfig.getParameter<std::vector<double> >("TS4TS5LowerThreshold");
0074   std::vector<double> TS4TS5LowerCutTemp = iConfig.getParameter<std::vector<double> >("TS4TS5LowerCut");
0075 
0076   for (int i = 0; i < (int)TS4TS5UpperThresholdTemp.size() && i < (int)TS4TS5UpperCutTemp.size(); i++)
0077     TS4TS5UpperCut_.push_back(std::pair<double, double>(TS4TS5UpperThresholdTemp[i], TS4TS5UpperCutTemp[i]));
0078   sort(TS4TS5UpperCut_.begin(), TS4TS5UpperCut_.end());
0079 
0080   for (int i = 0; i < (int)TS4TS5LowerThresholdTemp.size() && i < (int)TS4TS5LowerCutTemp.size(); i++)
0081     TS4TS5LowerCut_.push_back(std::pair<double, double>(TS4TS5LowerThresholdTemp[i], TS4TS5LowerCutTemp[i]));
0082   sort(TS4TS5LowerCut_.begin(), TS4TS5LowerCut_.end());
0083 
0084   m_theCaloMetToken = consumes<reco::CaloMETCollection>(CaloMetCollectionTag_);
0085   m_theHcalNoiseToken = consumes<reco::HcalNoiseRBXCollection>(HcalNoiseRBXCollectionTag_);
0086 
0087   produces<reco::CaloMETCollection>();
0088 }
0089 
0090 HLTHcalMETNoiseCleaner::~HLTHcalMETNoiseCleaner() = default;
0091 
0092 void HLTHcalMETNoiseCleaner::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0093   edm::ParameterSetDescription desc;
0094   desc.add<edm::InputTag>("HcalNoiseRBXCollection", edm::InputTag("hltHcalNoiseInfoProducer"));
0095   desc.add<edm::InputTag>("CaloMetCollection", edm::InputTag("hltMet"));
0096   desc.add<double>("CaloMetCut", 0.0);
0097   desc.add<int>("severity", 1);
0098   desc.add<int>("maxNumRBXs", 2);
0099   desc.add<int>("numRBXsToConsider", 2);
0100   desc.add<bool>("accept2NoiseRBXEvents", true);
0101   desc.add<bool>("needEMFCoincidence", true);
0102   desc.add<double>("minRBXEnergy", 50.0);
0103   desc.add<double>("minRatio", -999.);
0104   desc.add<double>("maxRatio", 999.);
0105   desc.add<int>("minHPDHits", 17);
0106   desc.add<int>("minRBXHits", 999);
0107   desc.add<int>("minHPDNoOtherHits", 10);
0108   desc.add<int>("minZeros", 10);
0109   desc.add<double>("minHighEHitTime", -9999.0);
0110   desc.add<double>("maxHighEHitTime", 9999.0);
0111   desc.add<double>("maxRBXEMF", 0.02);
0112   desc.add<double>("minRecHitE", 1.5);
0113   desc.add<double>("minLowHitE", 10.0);
0114   desc.add<double>("minHighHitE", 25.0);
0115   desc.add<double>("minR45HitE", 5.0);
0116   desc.add<double>("TS4TS5EnergyThreshold", 50.0);
0117 
0118   double TS4TS5UpperThresholdArray[5] = {70, 90, 100, 400, 4000};
0119   double TS4TS5UpperCutArray[5] = {1, 0.8, 0.75, 0.72, 0.72};
0120   double TS4TS5LowerThresholdArray[7] = {100, 120, 150, 200, 300, 400, 500};
0121   double TS4TS5LowerCutArray[7] = {-1, -0.7, -0.4, -0.2, -0.08, 0, 0.1};
0122   std::vector<double> TS4TS5UpperThreshold(TS4TS5UpperThresholdArray, TS4TS5UpperThresholdArray + 5);
0123   std::vector<double> TS4TS5UpperCut(TS4TS5UpperCutArray, TS4TS5UpperCutArray + 5);
0124   std::vector<double> TS4TS5LowerThreshold(TS4TS5LowerThresholdArray, TS4TS5LowerThresholdArray + 7);
0125   std::vector<double> TS4TS5LowerCut(TS4TS5LowerCutArray, TS4TS5LowerCutArray + 7);
0126 
0127   desc.add<std::vector<double> >("TS4TS5UpperThreshold", TS4TS5UpperThreshold);
0128   desc.add<std::vector<double> >("TS4TS5UpperCut", TS4TS5UpperCut);
0129   desc.add<std::vector<double> >("TS4TS5LowerThreshold", TS4TS5LowerThreshold);
0130   desc.add<std::vector<double> >("TS4TS5LowerCut", TS4TS5LowerCut);
0131   descriptions.add("hltHcalMETNoiseCleaner", desc);
0132 }
0133 
0134 //
0135 // member functions
0136 //
0137 
0138 bool HLTHcalMETNoiseCleaner::filter(edm::StreamID, edm::Event& iEvent, const edm::EventSetup& iSetup) const {
0139   using namespace reco;
0140 
0141   //output collection
0142   std::unique_ptr<CaloMETCollection> CleanedMET(new CaloMETCollection);
0143 
0144   //get the calo MET / MHT
0145   edm::Handle<CaloMETCollection> met_h;
0146   iEvent.getByToken(m_theCaloMetToken, met_h);
0147 
0148   if (not met_h.isValid() or met_h->empty() or
0149       met_h->front().pt() < 0) {  //No Valid MET, don't do anything and accept the event
0150     return true;                  // we shouldn't get here, but lets not crash
0151   }
0152 
0153   reco::CaloMET inCaloMet = met_h->front();
0154 
0155   // in this case, do not filter anything
0156   if (severity_ == 0) {
0157     CleanedMET->push_back(inCaloMet);
0158     iEvent.put(std::move(CleanedMET));
0159     return true;
0160   }
0161 
0162   // get the RBXs produced by RecoMET/METProducers/HcalNoiseInfoProducer
0163   edm::Handle<HcalNoiseRBXCollection> rbxs_h;
0164   iEvent.getByToken(m_theHcalNoiseToken, rbxs_h);
0165   if (!rbxs_h.isValid()) {
0166     edm::LogError("DataNotFound") << "HLTHcalMETNoiseCleaner: Could not find HcalNoiseRBXCollection product named "
0167                                   << HcalNoiseRBXCollectionTag_ << "." << std::endl;
0168     CleanedMET->push_back(inCaloMet);
0169     iEvent.put(std::move(CleanedMET));
0170     return true;  // no valid RBXs
0171   }
0172 
0173   // create a sorted set of the RBXs, ordered by energy
0174   noisedataset_t data;
0175   for (auto const& rbx : *rbxs_h) {
0176     CommonHcalNoiseRBXData d(rbx,
0177                              minRecHitE_,
0178                              minLowHitE_,
0179                              minHighHitE_,
0180                              TS4TS5EnergyThreshold_,
0181                              TS4TS5UpperCut_,
0182                              TS4TS5LowerCut_,
0183                              minR45HitE_);
0184     data.insert(d);
0185   }
0186   //if 0 RBXs are in the list, just accept
0187   if (data.empty()) {
0188     CleanedMET->push_back(inCaloMet);
0189     iEvent.put(std::move(CleanedMET));
0190     return true;
0191   }
0192   // data is now sorted by RBX energy
0193   // only consider top N=numRBXsToConsider_ energy RBXs
0194   int cntr = 0;
0195   int nNoise = 0;
0196 
0197   TVector3 metVec;
0198   metVec.SetPtEtaPhi(met_h->front().pt(), 0, met_h->front().phi());
0199 
0200   TVector3 noiseHPDVector(0, 0, 0);
0201   TVector3 secondHPDVector(0, 0, 0);
0202   for (auto it = data.begin(); it != data.end() && cntr < numRBXsToConsider_; it++, cntr++) {
0203     bool isNoise = false;
0204     bool passFilter = true;
0205     bool passEMF = true;
0206     if (it->energy() > minRBXEnergy_) {
0207       if (it->validRatio() && it->ratio() < minRatio_)
0208         passFilter = false;
0209       else if (it->validRatio() && it->ratio() > maxRatio_)
0210         passFilter = false;
0211       else if (it->numHPDHits() >= minHPDHits_)
0212         passFilter = false;
0213       else if (it->numRBXHits() >= minRBXHits_)
0214         passFilter = false;
0215       else if (it->numHPDNoOtherHits() >= minHPDNoOtherHits_)
0216         passFilter = false;
0217       else if (it->numZeros() >= minZeros_)
0218         passFilter = false;
0219       else if (it->minHighEHitTime() < minHighEHitTime_)
0220         passFilter = false;
0221       else if (it->maxHighEHitTime() > maxHighEHitTime_)
0222         passFilter = false;
0223       else if (!it->PassTS4TS5())
0224         passFilter = false;
0225 
0226       if (it->RBXEMF() < maxRBXEMF_) {
0227         passEMF = false;
0228       }
0229     }
0230 
0231     if ((needEMFCoincidence_ && !passEMF && !passFilter) || (!needEMFCoincidence_ && !passFilter)) {  // check for noise
0232       LogDebug("") << "HLTHcalMETNoiseCleaner debug: Found a noisy RBX: "
0233                    << "energy=" << it->energy() << "; "
0234                    << "ratio=" << it->ratio() << "; "
0235                    << "# RBX hits=" << it->numRBXHits() << "; "
0236                    << "# HPD hits=" << it->numHPDHits() << "; "
0237                    << "# Zeros=" << it->numZeros() << "; "
0238                    << "min time=" << it->minHighEHitTime() << "; "
0239                    << "max time=" << it->maxHighEHitTime() << "; "
0240                    << "passTS4TS5=" << it->PassTS4TS5() << "; "
0241                    << "RBX EMF=" << it->RBXEMF() << std::endl;
0242       nNoise++;
0243       isNoise = true;
0244     }  // OK, checked for noise
0245 
0246     //------------First Noisy RBX-----------------------
0247     if (isNoise && nNoise == 1) {
0248       edm::RefVector<CaloTowerCollection> noiseTowers = it->rbxTowers();
0249       edm::RefVector<CaloTowerCollection>::const_iterator noiseTowersIt;
0250       // get the energy vector for this RBX from the calotowers
0251       for (noiseTowersIt = noiseTowers.begin(); noiseTowersIt != noiseTowers.end(); noiseTowersIt++) {
0252         TVector3 towerVec;
0253         towerVec.SetPtEtaPhi((*noiseTowersIt)->pt(), (*noiseTowersIt)->eta(), (*noiseTowersIt)->phi());
0254         noiseHPDVector += towerVec;  // add this tower to the vector for the RBX
0255       }
0256       if (noiseHPDVector.Mag() > 0)
0257         noiseHPDVector.SetPtEtaPhi(noiseHPDVector.Pt(), 0, noiseHPDVector.Phi());  // make the noise transverse
0258       else
0259         noiseHPDVector.SetPtEtaPhi(0, 0, 0);
0260     }
0261     //-----------FOUND a SECOND NOISY RBX-------------------
0262     if (isNoise && cntr > 0) {
0263       CleanedMET->push_back(inCaloMet);
0264       iEvent.put(std::move(CleanedMET));
0265       return accept2NoiseRBXEvents_;  // don't try to clean these for the moment, just keep or throw away
0266     }
0267     //----------LEADING RBX is NOT NOISY--------------------
0268     if (!isNoise && cntr == 0) {
0269       CleanedMET->push_back(inCaloMet);
0270       iEvent.put(std::move(CleanedMET));
0271       return true;  // don't reject the event if the leading RBX isn't noise
0272     }
0273     //-----------SUBLEADING RBX is NOT NOISY: STORE INFO----
0274     if (!isNoise && nNoise > 0) {  //second RBX isn't noisy (and first one was), so clean
0275       edm::RefVector<CaloTowerCollection> noiseTowers = it->rbxTowers();
0276       edm::RefVector<CaloTowerCollection>::const_iterator noiseTowersIt;
0277       for (noiseTowersIt = noiseTowers.begin(); noiseTowersIt != noiseTowers.end(); noiseTowersIt++) {
0278         TVector3 towerVec;
0279         towerVec.SetPtEtaPhi((*noiseTowersIt)->pt(), (*noiseTowersIt)->eta(), (*noiseTowersIt)->phi());
0280         secondHPDVector += towerVec;  // add this tower to the vector for the RBX
0281       }
0282       if (secondHPDVector.Mag() > 0)
0283         secondHPDVector.SetPtEtaPhi(secondHPDVector.Pt(), 0, secondHPDVector.Phi());  // make the second transverse
0284       else
0285         secondHPDVector.SetPtEtaPhi(0, 0, 0);
0286       break;
0287     }
0288   }  // end RBX loop
0289 
0290   if (noiseHPDVector.Mag() == 0) {
0291     CleanedMET->push_back(inCaloMet);
0292     iEvent.put(std::move(CleanedMET));
0293     return true;  // don't reject the event if the leading RBX isn't noise
0294   }
0295 
0296   //********************************************************************************
0297   //The Event gets here only if it had exactly 1 noisy RBX in the lead position
0298   //********************************************************************************
0299 
0300   float METsumet = met_h->front().energy();
0301 
0302   metVec += noiseHPDVector;
0303 
0304   float ZMETsumet = METsumet - noiseHPDVector.Mag();
0305   float ZMETpt = metVec.Pt();
0306   float ZMETphi = metVec.Phi();
0307 
0308   //put the second RBX vector in the eta phi position of the leading RBX vector
0309 
0310   float SMETsumet = 0;
0311   float SMETpt = 0;
0312   float SMETphi = 0;
0313   if (secondHPDVector.Mag() > 0.) {
0314     secondHPDVector.SetPtEtaPhi(secondHPDVector.Pt(), noiseHPDVector.Eta(), noiseHPDVector.Phi());
0315     metVec -= secondHPDVector;
0316     SMETsumet = METsumet - noiseHPDVector.Mag();
0317     SMETpt = metVec.Pt();
0318     SMETphi = metVec.Phi();
0319   }
0320   //Get the maximum MET:
0321   float CorMetSumEt, CorMetPt, CorMetPhi;
0322   if (ZMETpt > SMETpt) {
0323     CorMetSumEt = ZMETsumet;
0324     CorMetPt = ZMETpt;
0325     CorMetPhi = ZMETphi;
0326   } else {
0327     CorMetSumEt = SMETsumet;
0328     CorMetPt = SMETpt;
0329     CorMetPhi = SMETphi;
0330   }
0331 
0332   reco::CaloMET corMet = BuildCaloMet(CorMetSumEt, CorMetPt, CorMetPhi);
0333   CleanedMET->push_back(corMet);
0334   iEvent.put(std::move(CleanedMET));
0335 
0336   return (corMet.pt() > CaloMetCut_);
0337 }
0338 
0339 reco::CaloMET HLTHcalMETNoiseCleaner::BuildCaloMet(float sumet, float pt, float phi) const {
0340   // Instantiate the container to hold the calorimeter specific information
0341 
0342   typedef math::XYZPoint Point;
0343   typedef math::XYZTLorentzVector LorentzVector;
0344 
0345   SpecificCaloMETData specific;
0346   // Initialise the container
0347   specific.MaxEtInEmTowers = 0.0;     // Maximum energy in EM towers
0348   specific.MaxEtInHadTowers = 0.0;    // Maximum energy in HCAL towers
0349   specific.HadEtInHO = 0.0;           // Hadronic energy fraction in HO
0350   specific.HadEtInHB = 0.0;           // Hadronic energy in HB
0351   specific.HadEtInHF = 0.0;           // Hadronic energy in HF
0352   specific.HadEtInHE = 0.0;           // Hadronic energy in HE
0353   specific.EmEtInEB = 0.0;            // Em energy in EB
0354   specific.EmEtInEE = 0.0;            // Em energy in EE
0355   specific.EmEtInHF = 0.0;            // Em energy in HF
0356   specific.EtFractionHadronic = 0.0;  // Hadronic energy fraction
0357   specific.EtFractionEm = 0.0;        // Em energy fraction
0358   specific.CaloSETInpHF = 0.0;        // CaloSET in HF+
0359   specific.CaloSETInmHF = 0.0;        // CaloSET in HF-
0360   specific.CaloMETInpHF = 0.0;        // CaloMET in HF+
0361   specific.CaloMETInmHF = 0.0;        // CaloMET in HF-
0362   specific.CaloMETPhiInpHF = 0.0;     // CaloMET-phi in HF+
0363   specific.CaloMETPhiInmHF = 0.0;     // CaloMET-phi in HF-
0364   specific.METSignificance = 0.0;
0365 
0366   TLorentzVector p4TL;
0367   p4TL.SetPtEtaPhiM(pt, 0., phi, 0.);
0368   const LorentzVector p4(p4TL.X(), p4TL.Y(), 0, p4TL.T());
0369   const Point vtx(0.0, 0.0, 0.0);
0370   reco::CaloMET specificmet(specific, sumet, p4, vtx);
0371   return specificmet;
0372 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!