File indexing completed on 2024-04-06 12:18:33
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
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
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
0136
0137
0138 bool HLTHcalMETNoiseCleaner::filter(edm::StreamID, edm::Event& iEvent, const edm::EventSetup& iSetup) const {
0139 using namespace reco;
0140
0141
0142 std::unique_ptr<CaloMETCollection> CleanedMET(new CaloMETCollection);
0143
0144
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) {
0150 return true;
0151 }
0152
0153 reco::CaloMET inCaloMet = met_h->front();
0154
0155
0156 if (severity_ == 0) {
0157 CleanedMET->push_back(inCaloMet);
0158 iEvent.put(std::move(CleanedMET));
0159 return true;
0160 }
0161
0162
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;
0171 }
0172
0173
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
0187 if (data.empty()) {
0188 CleanedMET->push_back(inCaloMet);
0189 iEvent.put(std::move(CleanedMET));
0190 return true;
0191 }
0192
0193
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)) {
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 }
0245
0246
0247 if (isNoise && nNoise == 1) {
0248 edm::RefVector<CaloTowerCollection> noiseTowers = it->rbxTowers();
0249 edm::RefVector<CaloTowerCollection>::const_iterator noiseTowersIt;
0250
0251 for (noiseTowersIt = noiseTowers.begin(); noiseTowersIt != noiseTowers.end(); noiseTowersIt++) {
0252 TVector3 towerVec;
0253 towerVec.SetPtEtaPhi((*noiseTowersIt)->pt(), (*noiseTowersIt)->eta(), (*noiseTowersIt)->phi());
0254 noiseHPDVector += towerVec;
0255 }
0256 if (noiseHPDVector.Mag() > 0)
0257 noiseHPDVector.SetPtEtaPhi(noiseHPDVector.Pt(), 0, noiseHPDVector.Phi());
0258 else
0259 noiseHPDVector.SetPtEtaPhi(0, 0, 0);
0260 }
0261
0262 if (isNoise && cntr > 0) {
0263 CleanedMET->push_back(inCaloMet);
0264 iEvent.put(std::move(CleanedMET));
0265 return accept2NoiseRBXEvents_;
0266 }
0267
0268 if (!isNoise && cntr == 0) {
0269 CleanedMET->push_back(inCaloMet);
0270 iEvent.put(std::move(CleanedMET));
0271 return true;
0272 }
0273
0274 if (!isNoise && nNoise > 0) {
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;
0281 }
0282 if (secondHPDVector.Mag() > 0)
0283 secondHPDVector.SetPtEtaPhi(secondHPDVector.Pt(), 0, secondHPDVector.Phi());
0284 else
0285 secondHPDVector.SetPtEtaPhi(0, 0, 0);
0286 break;
0287 }
0288 }
0289
0290 if (noiseHPDVector.Mag() == 0) {
0291 CleanedMET->push_back(inCaloMet);
0292 iEvent.put(std::move(CleanedMET));
0293 return true;
0294 }
0295
0296
0297
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
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
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
0341
0342 typedef math::XYZPoint Point;
0343 typedef math::XYZTLorentzVector LorentzVector;
0344
0345 SpecificCaloMETData specific;
0346
0347 specific.MaxEtInEmTowers = 0.0;
0348 specific.MaxEtInHadTowers = 0.0;
0349 specific.HadEtInHO = 0.0;
0350 specific.HadEtInHB = 0.0;
0351 specific.HadEtInHF = 0.0;
0352 specific.HadEtInHE = 0.0;
0353 specific.EmEtInEB = 0.0;
0354 specific.EmEtInEE = 0.0;
0355 specific.EmEtInHF = 0.0;
0356 specific.EtFractionHadronic = 0.0;
0357 specific.EtFractionEm = 0.0;
0358 specific.CaloSETInpHF = 0.0;
0359 specific.CaloSETInmHF = 0.0;
0360 specific.CaloMETInpHF = 0.0;
0361 specific.CaloMETInmHF = 0.0;
0362 specific.CaloMETPhiInpHF = 0.0;
0363 specific.CaloMETPhiInmHF = 0.0;
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 }