Back to home page

Project CMSSW displayed by LXR
 
 

    

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

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