File indexing completed on 2021-02-14 14:20:31
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_(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
0139
0140
0141 bool HLTHcalMETNoiseCleaner::filter(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0142 using namespace reco;
0143
0144
0145 std::unique_ptr<CaloMETCollection> CleanedMET(new CaloMETCollection);
0146
0147
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) {
0153 return true;
0154 }
0155
0156 reco::CaloMET inCaloMet = met_h->front();
0157
0158
0159 if (severity_ == 0) {
0160 CleanedMET->push_back(inCaloMet);
0161 iEvent.put(std::move(CleanedMET));
0162 return true;
0163 }
0164
0165
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;
0174 }
0175
0176
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
0190 if (data.empty()) {
0191 CleanedMET->push_back(inCaloMet);
0192 iEvent.put(std::move(CleanedMET));
0193 return true;
0194 }
0195
0196
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)) {
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 }
0248
0249
0250 if (isNoise && nNoise == 1) {
0251 edm::RefVector<CaloTowerCollection> noiseTowers = it->rbxTowers();
0252 edm::RefVector<CaloTowerCollection>::const_iterator noiseTowersIt;
0253
0254 for (noiseTowersIt = noiseTowers.begin(); noiseTowersIt != noiseTowers.end(); noiseTowersIt++) {
0255 TVector3 towerVec;
0256 towerVec.SetPtEtaPhi((*noiseTowersIt)->pt(), (*noiseTowersIt)->eta(), (*noiseTowersIt)->phi());
0257 noiseHPDVector += towerVec;
0258 }
0259 if (noiseHPDVector.Mag() > 0)
0260 noiseHPDVector.SetPtEtaPhi(noiseHPDVector.Pt(), 0, noiseHPDVector.Phi());
0261 else
0262 noiseHPDVector.SetPtEtaPhi(0, 0, 0);
0263 }
0264
0265 if (isNoise && cntr > 0) {
0266 CleanedMET->push_back(inCaloMet);
0267 iEvent.put(std::move(CleanedMET));
0268 return accept2NoiseRBXEvents_;
0269 }
0270
0271 if (!isNoise && cntr == 0) {
0272 CleanedMET->push_back(inCaloMet);
0273 iEvent.put(std::move(CleanedMET));
0274 return true;
0275 }
0276
0277 if (!isNoise && nNoise > 0) {
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;
0284 }
0285 if (secondHPDVector.Mag() > 0)
0286 secondHPDVector.SetPtEtaPhi(secondHPDVector.Pt(), 0, secondHPDVector.Phi());
0287 else
0288 secondHPDVector.SetPtEtaPhi(0, 0, 0);
0289 break;
0290 }
0291 }
0292
0293 if (noiseHPDVector.Mag() == 0) {
0294 CleanedMET->push_back(inCaloMet);
0295 iEvent.put(std::move(CleanedMET));
0296 return true;
0297 }
0298
0299
0300
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
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
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
0344
0345 typedef math::XYZPoint Point;
0346 typedef math::XYZTLorentzVector LorentzVector;
0347
0348 SpecificCaloMETData specific;
0349
0350 specific.MaxEtInEmTowers = 0.0;
0351 specific.MaxEtInHadTowers = 0.0;
0352 specific.HadEtInHO = 0.0;
0353 specific.HadEtInHB = 0.0;
0354 specific.HadEtInHF = 0.0;
0355 specific.HadEtInHE = 0.0;
0356 specific.EmEtInEB = 0.0;
0357 specific.EmEtInEE = 0.0;
0358 specific.EmEtInHF = 0.0;
0359 specific.EtFractionHadronic = 0.0;
0360 specific.EtFractionEm = 0.0;
0361 specific.CaloSETInpHF = 0.0;
0362 specific.CaloSETInmHF = 0.0;
0363 specific.CaloMETInpHF = 0.0;
0364 specific.CaloMETInmHF = 0.0;
0365 specific.CaloMETPhiInpHF = 0.0;
0366 specific.CaloMETPhiInmHF = 0.0;
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 }