diff --git a/DataFormats/METReco/interface/HcalHaloData.h b/DataFormats/METReco/interface/HcalHaloData.h index 93e92e68afd81..3cec018dd9d52 100644 --- a/DataFormats/METReco/interface/HcalHaloData.h +++ b/DataFormats/METReco/interface/HcalHaloData.h @@ -14,13 +14,19 @@ struct HaloTowerStrip { std::vector > cellTowerIds; float hadEt; + float energyRatio; + float emEt; HaloTowerStrip() { hadEt = 0.0; + energyRatio = -1.0; + emEt = 0.0; cellTowerIds.clear(); } HaloTowerStrip(const HaloTowerStrip& strip) { cellTowerIds = strip.cellTowerIds; hadEt = strip.hadEt; + energyRatio = strip.energyRatio; + emEt = strip.emEt; } }; diff --git a/DataFormats/METReco/src/classes_def.xml b/DataFormats/METReco/src/classes_def.xml index c339d2dd3614d..5a3ac40872017 100644 --- a/DataFormats/METReco/src/classes_def.xml +++ b/DataFormats/METReco/src/classes_def.xml @@ -172,7 +172,8 @@ - + + diff --git a/RecoMET/METAlgorithms/src/HcalHaloAlgo.cc b/RecoMET/METAlgorithms/src/HcalHaloAlgo.cc index 224d50f366101..efea2bd7ca639 100644 --- a/RecoMET/METAlgorithms/src/HcalHaloAlgo.cc +++ b/RecoMET/METAlgorithms/src/HcalHaloAlgo.cc @@ -1,4 +1,5 @@ #include "RecoMET/METAlgorithms/interface/HcalHaloAlgo.h" +#include /* [class]: HcalHaloAlgo @@ -128,24 +129,38 @@ HcalHaloData HcalHaloAlgo::Calculate(const CaloGeometry& TheCaloGeometry, edm::H } } + // Don't use HF. int maxAbsIEta = 29; + + + std::map iPhiHadEtMap; std::vector sortedCaloTowers; for(CaloTowerCollection::const_iterator tower = TheCaloTowers->begin(); tower != TheCaloTowers->end(); tower++) { - if(abs(tower->ieta()) <= maxAbsIEta && tower->numProblematicHcalCells() > 0) - sortedCaloTowers.push_back(&(*tower)); + if(abs(tower->ieta()) > maxAbsIEta) continue; + + int iPhi = tower->iphi(); + if(!iPhiHadEtMap.count(iPhi)) iPhiHadEtMap[iPhi] = 0.0; + iPhiHadEtMap[iPhi] += tower->hadEt(); + + if(tower->numProblematicHcalCells() > 0) sortedCaloTowers.push_back(&(*tower)); + } + // Sort towers such that lowest iphi and ieta are first, highest last, and towers // with same iphi value are consecutive. Then we can do everything else in one loop. std::sort(sortedCaloTowers.begin(), sortedCaloTowers.end(), CompareTowers); HaloTowerStrip strip; + int prevIEta = -99, prevIPhi = -99; float prevHadEt = 0.; + float prevEmEt = 0.; std::pair prevPair, towerPair; bool wasContiguous = true; + // Loop through and store a vector of pairs (problematicCells, DetId) for each contiguous strip we find for(unsigned int i = 0; i < sortedCaloTowers.size(); i++) { const CaloTower* tower = sortedCaloTowers[i]; @@ -162,22 +177,40 @@ HcalHaloData HcalHaloAlgo::Calculate(const CaloGeometry& TheCaloGeometry, edm::H strip.cellTowerIds.push_back(prevPair); strip.cellTowerIds.push_back(towerPair); strip.hadEt += prevHadEt + tower->hadEt(); + strip.emEt += prevEmEt + tower->emEt(); } if(wasContiguous && isContiguous) { strip.cellTowerIds.push_back(towerPair); strip.hadEt += tower->hadEt(); + strip.emEt += tower->emEt(); } if((wasContiguous && !isContiguous) || i == sortedCaloTowers.size()-1) { //ended the strip, so flush it - if(strip.cellTowerIds.size() > 2) { + + if(strip.cellTowerIds.size() > 3) { + + int iPhi = strip.cellTowerIds.at(0).second.iphi(); + int iPhiLower = (iPhi == 1) ? 72 : iPhi - 1; + int iPhiUpper = (iPhi == 72) ? 1 : iPhi + 1; + + float energyRatio = 0.0; + if(iPhiHadEtMap.count(iPhiLower)) energyRatio += iPhiHadEtMap[iPhiLower]; + if(iPhiHadEtMap.count(iPhiUpper)) energyRatio += iPhiHadEtMap[iPhiUpper]; + iPhiHadEtMap[iPhi] = max(iPhiHadEtMap[iPhi], 0.001F); + + energyRatio /= iPhiHadEtMap[iPhi]; + strip.energyRatio = energyRatio; + TheHcalHaloData.getProblematicStrips().push_back( strip ); + } strip = HaloTowerStrip(); } wasContiguous = isContiguous; prevPair = towerPair; + prevEmEt = tower->emEt(); prevIPhi = tower->iphi(); prevIEta = tower->ieta(); prevHadEt = tower->hadEt(); diff --git a/RecoMET/METFilters/plugins/HcalStripHaloFilter.cc b/RecoMET/METFilters/plugins/HcalStripHaloFilter.cc index a3d4e0a88f08f..ccd3718fb973e 100644 --- a/RecoMET/METFilters/plugins/HcalStripHaloFilter.cc +++ b/RecoMET/METFilters/plugins/HcalStripHaloFilter.cc @@ -18,12 +18,16 @@ class HcalStripHaloFilter : public edm::global::EDFilter<> { const bool taggingMode_; const int maxWeightedStripLength_; + const double maxEnergyRatio_; + const double minHadEt_; edm::EDGetTokenT beamHaloSummaryToken_; }; HcalStripHaloFilter::HcalStripHaloFilter(const edm::ParameterSet & iConfig) : taggingMode_ (iConfig.getParameter ("taggingMode")) , maxWeightedStripLength_ (iConfig.getParameter("maxWeightedStripLength")) + , maxEnergyRatio_ (iConfig.getParameter("maxEnergyRatio")) + , minHadEt_ (iConfig.getParameter("minHadEt")) , beamHaloSummaryToken_(consumes(edm::InputTag("BeamHaloSummary"))) { produces(); @@ -42,7 +46,9 @@ bool HcalStripHaloFilter::filter(edm::StreamID iID, edm::Event & iEvent, const e for (unsigned int iTower = 0; iTower < problematicStrip.cellTowerIds.size(); iTower++) { numContiguousCells += (int)problematicStrip.cellTowerIds[iTower].first; } - if(numContiguousCells > maxWeightedStripLength_) { + if( numContiguousCells > maxWeightedStripLength_ + && problematicStrip.energyRatio < maxEnergyRatio_ + && problematicStrip.hadEt > minHadEt_ ) { pass = false; break; } diff --git a/RecoMET/METFilters/python/HcalStripHaloFilter_cfi.py b/RecoMET/METFilters/python/HcalStripHaloFilter_cfi.py index a297d79160795..f89a2b5f48612 100644 --- a/RecoMET/METFilters/python/HcalStripHaloFilter_cfi.py +++ b/RecoMET/METFilters/python/HcalStripHaloFilter_cfi.py @@ -3,5 +3,7 @@ HcalStripHaloFilter = cms.EDFilter( "HcalStripHaloFilter", taggingMode = cms.bool(False), - maxWeightedStripLength = cms.int32(9) # values higher than this are rejected + maxWeightedStripLength = cms.int32(7), + maxEnergyRatio = cms.double(0.15), + minHadEt = cms.double(100.0) )