File indexing completed on 2022-01-04 23:06:06
0001 #include "HLTEcalResonanceFilter.h"
0002 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0003 #include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
0004 #include "TLorentzVector.h"
0005
0006 using namespace std;
0007 using namespace edm;
0008
0009 HLTEcalResonanceFilter::HLTEcalResonanceFilter(const edm::ParameterSet &iConfig)
0010 : caloTopologyRecordToken_(esConsumes()),
0011 ecalChannelStatusRcdToken_(esConsumes()),
0012 caloGeometryRecordToken_(esConsumes()) {
0013 barrelHits_ = iConfig.getParameter<edm::InputTag>("barrelHits");
0014 barrelClusters_ = iConfig.getParameter<edm::InputTag>("barrelClusters");
0015 barrelHitsToken_ = consumes<EBRecHitCollection>(barrelHits_);
0016 barrelClustersToken_ = consumes<reco::BasicClusterCollection>(barrelClusters_);
0017
0018 endcapHits_ = iConfig.getParameter<edm::InputTag>("endcapHits");
0019 endcapClusters_ = iConfig.getParameter<edm::InputTag>("endcapClusters");
0020 endcapHitsToken_ = consumes<EERecHitCollection>(endcapHits_);
0021 endcapClustersToken_ = consumes<reco::BasicClusterCollection>(endcapClusters_);
0022
0023 doSelBarrel_ = iConfig.getParameter<bool>("doSelBarrel");
0024 if (doSelBarrel_) {
0025 edm::ParameterSet barrelSelection = iConfig.getParameter<edm::ParameterSet>("barrelSelection");
0026
0027
0028 selePtGamma_ = barrelSelection.getParameter<double>("selePtGamma");
0029 selePtPair_ = barrelSelection.getParameter<double>("selePtPair");
0030 seleS4S9Gamma_ = barrelSelection.getParameter<double>("seleS4S9Gamma");
0031 seleS9S25Gamma_ = barrelSelection.getParameter<double>("seleS9S25Gamma");
0032 seleMinvMaxBarrel_ = barrelSelection.getParameter<double>("seleMinvMaxBarrel");
0033 seleMinvMinBarrel_ = barrelSelection.getParameter<double>("seleMinvMinBarrel");
0034 ptMinForIsolation_ = barrelSelection.getParameter<double>("ptMinForIsolation");
0035 removePi0CandidatesForEta_ = barrelSelection.getParameter<bool>("removePi0CandidatesForEta");
0036 if (removePi0CandidatesForEta_) {
0037 massLowPi0Cand_ = barrelSelection.getParameter<double>("massLowPi0Cand");
0038 massHighPi0Cand_ = barrelSelection.getParameter<double>("massHighPi0Cand");
0039 }
0040 seleIso_ = barrelSelection.getParameter<double>("seleIso");
0041 ptMinForIsolation_ = barrelSelection.getParameter<double>("ptMinForIsolation");
0042 seleBeltDR_ = barrelSelection.getParameter<double>("seleBeltDR");
0043 seleBeltDeta_ = barrelSelection.getParameter<double>("seleBeltDeta");
0044
0045 store5x5RecHitEB_ = barrelSelection.getParameter<bool>("store5x5RecHitEB");
0046 BarrelHits_ = barrelSelection.getParameter<string>("barrelHitCollection");
0047
0048 produces<EBRecHitCollection>(BarrelHits_);
0049 }
0050
0051 doSelEndcap_ = iConfig.getParameter<bool>("doSelEndcap");
0052 if (doSelEndcap_) {
0053 edm::ParameterSet endcapSelection = iConfig.getParameter<edm::ParameterSet>("endcapSelection");
0054
0055
0056 seleMinvMaxEndCap_ = endcapSelection.getParameter<double>("seleMinvMaxEndCap");
0057 seleMinvMinEndCap_ = endcapSelection.getParameter<double>("seleMinvMinEndCap");
0058
0059 region1_EndCap_ = endcapSelection.getParameter<double>("region1_EndCap");
0060 selePtGammaEndCap_region1_ = endcapSelection.getParameter<double>("selePtGammaEndCap_region1");
0061 selePtPairEndCap_region1_ = endcapSelection.getParameter<double>("selePtPairEndCap_region1");
0062
0063 region2_EndCap_ = endcapSelection.getParameter<double>("region2_EndCap");
0064 selePtGammaEndCap_region2_ = endcapSelection.getParameter<double>("selePtGammaEndCap_region2");
0065 selePtPairEndCap_region2_ = endcapSelection.getParameter<double>("selePtPairEndCap_region2");
0066
0067 selePtGammaEndCap_region3_ = endcapSelection.getParameter<double>("selePtGammaEndCap_region3");
0068 selePtPairEndCap_region3_ = endcapSelection.getParameter<double>("selePtPairEndCap_region3");
0069 selePtPairMaxEndCap_region3_ = endcapSelection.getParameter<double>("selePtPairMaxEndCap_region3");
0070
0071 seleS4S9GammaEndCap_ = endcapSelection.getParameter<double>("seleS4S9GammaEndCap");
0072 seleS9S25GammaEndCap_ = endcapSelection.getParameter<double>("seleS9S25GammaEndCap");
0073 ptMinForIsolationEndCap_ = endcapSelection.getParameter<double>("ptMinForIsolationEndCap");
0074 seleBeltDREndCap_ = endcapSelection.getParameter<double>("seleBeltDREndCap");
0075 seleBeltDetaEndCap_ = endcapSelection.getParameter<double>("seleBeltDetaEndCap");
0076 seleIsoEndCap_ = endcapSelection.getParameter<double>("seleIsoEndCap");
0077
0078 store5x5RecHitEE_ = endcapSelection.getParameter<bool>("store5x5RecHitEE");
0079 EndcapHits_ = endcapSelection.getParameter<string>("endcapHitCollection");
0080 produces<EERecHitCollection>(EndcapHits_);
0081 }
0082
0083 useRecoFlag_ = iConfig.getParameter<bool>("useRecoFlag");
0084 flagLevelRecHitsToUse_ = iConfig.getParameter<int>("flagLevelRecHitsToUse");
0085
0086 useDBStatus_ = iConfig.getParameter<bool>("useDBStatus");
0087 statusLevelRecHitsToUse_ = iConfig.getParameter<int>("statusLevelRecHitsToUse");
0088
0089 preshHitProducer_ = iConfig.getParameter<edm::InputTag>("preshRecHitProducer");
0090 preshHitsToken_ = consumes<EBRecHitCollection>(preshHitProducer_);
0091
0092
0093 storeRecHitES_ = iConfig.getParameter<bool>("storeRecHitES");
0094 if (storeRecHitES_) {
0095 edm::ParameterSet preshowerSelection = iConfig.getParameter<edm::ParameterSet>("preshowerSelection");
0096
0097
0098 preshNclust_ = preshowerSelection.getParameter<int>("preshNclust");
0099
0100 preshClustECut = preshowerSelection.getParameter<double>("preshClusterEnergyCut");
0101
0102 float preshStripECut = preshowerSelection.getParameter<double>("preshStripEnergyCut");
0103 int preshSeededNst = preshowerSelection.getParameter<int>("preshSeededNstrip");
0104
0105 calib_planeX_ = preshowerSelection.getParameter<double>("preshCalibPlaneX");
0106 calib_planeY_ = preshowerSelection.getParameter<double>("preshCalibPlaneY");
0107 gamma_ = preshowerSelection.getParameter<double>("preshCalibGamma");
0108 mip_ = preshowerSelection.getParameter<double>("preshCalibMIP");
0109
0110
0111 presh_algo_ = std::make_unique<PreshowerClusterAlgo>(preshStripECut, preshClustECut, preshSeededNst);
0112
0113 ESHits_ = preshowerSelection.getParameter<std::string>("ESCollection");
0114 produces<ESRecHitCollection>(ESHits_);
0115 }
0116
0117 debug_ = iConfig.getParameter<int>("debugLevel");
0118 }
0119
0120 HLTEcalResonanceFilter::~HLTEcalResonanceFilter() {}
0121
0122 void HLTEcalResonanceFilter::fillDescriptions(edm::ConfigurationDescriptions &descriptions) {
0123 edm::ParameterSetDescription desc;
0124 desc.add<edm::InputTag>("barrelHits", edm::InputTag("hltEcalRegionalPi0EtaRecHit", "EcalRecHitsEB"));
0125 desc.add<edm::InputTag>("endcapHits", edm::InputTag("hltEcalRegionalPi0EtaRecHit", "EcalRecHitsEE"));
0126 desc.add<edm::InputTag>("preshRecHitProducer", edm::InputTag("hltEcalRegionalPi0EtaRecHit", "EcalRecHitsES"));
0127 desc.add<edm::InputTag>("barrelClusters", edm::InputTag("hltSimple3x3Clusters", "Simple3x3ClustersBarrel"));
0128 desc.add<edm::InputTag>("endcapClusters", edm::InputTag("hltSimple3x3Clusters", "Simple3x3ClustersEndcap"));
0129 desc.add<bool>("useRecoFlag", false);
0130 desc.add<int>("flagLevelRecHitsToUse", 1);
0131 desc.add<bool>("useDBStatus", true);
0132 desc.add<int>("statusLevelRecHitsToUse", 1);
0133
0134 desc.add<bool>("doSelBarrel", true);
0135 edm::ParameterSetDescription barrelSelection;
0136 barrelSelection.add<double>("selePtGamma", 1.);
0137 barrelSelection.add<double>("selePtPair", 2.);
0138 barrelSelection.add<double>("seleMinvMaxBarrel", 0.22);
0139 barrelSelection.add<double>("seleMinvMinBarrel", 0.06);
0140 barrelSelection.add<bool>("removePi0CandidatesForEta", false);
0141 barrelSelection.add<double>("massLowPi0Cand", 0.104);
0142 barrelSelection.add<double>("massHighPi0Cand", 0.163);
0143 barrelSelection.add<double>("seleS4S9Gamma", 0.83);
0144 barrelSelection.add<double>("seleS9S25Gamma", 0.);
0145 barrelSelection.add<double>("ptMinForIsolation", 1.0);
0146 barrelSelection.add<double>("seleIso", 0.5);
0147 barrelSelection.add<double>("seleBeltDR", 0.2);
0148 barrelSelection.add<double>("seleBeltDeta", 0.05);
0149 barrelSelection.add<bool>("store5x5RecHitEB", false);
0150 barrelSelection.add<std::string>("barrelHitCollection", "pi0EcalRecHitsEB");
0151 desc.add<edm::ParameterSetDescription>("barrelSelection", barrelSelection);
0152
0153 desc.add<bool>("doSelEndcap", true);
0154 edm::ParameterSetDescription endcapSelection;
0155 endcapSelection.add<double>("seleMinvMaxEndCap", 0.3);
0156 endcapSelection.add<double>("seleMinvMinEndCap", 0.05);
0157 endcapSelection.add<double>("region1_EndCap", 2.0);
0158 endcapSelection.add<double>("selePtGammaEndCap_region1", 0.8);
0159 endcapSelection.add<double>("selePtPairEndCap_region1", 3.0);
0160 endcapSelection.add<double>("region2_EndCap", 2.5);
0161 endcapSelection.add<double>("selePtGammaEndCap_region2", 0.5);
0162 endcapSelection.add<double>("selePtPairEndCap_region2", 2.0);
0163 endcapSelection.add<double>("selePtGammaEndCap_region3", 0.3);
0164 endcapSelection.add<double>("selePtPairEndCap_region3", 1.2);
0165 endcapSelection.add<double>("selePtPairMaxEndCap_region3", 2.5);
0166 endcapSelection.add<double>("seleS4S9GammaEndCap", 0.9);
0167 endcapSelection.add<double>("seleS9S25GammaEndCap", 0.);
0168 endcapSelection.add<double>("ptMinForIsolationEndCap", 0.5);
0169 endcapSelection.add<double>("seleIsoEndCap", 0.5);
0170 endcapSelection.add<double>("seleBeltDREndCap", 0.2);
0171 endcapSelection.add<double>("seleBeltDetaEndCap", 0.05);
0172 endcapSelection.add<bool>("store5x5RecHitEE", false);
0173 endcapSelection.add<std::string>("endcapHitCollection", "pi0EcalRecHitsEE");
0174 desc.add<edm::ParameterSetDescription>("endcapSelection", endcapSelection);
0175
0176 desc.add<bool>("storeRecHitES", true);
0177 edm::ParameterSetDescription preshowerSelection;
0178 preshowerSelection.add<std::string>("ESCollection", "pi0EcalRecHitsES");
0179 preshowerSelection.add<int>("preshNclust", 4);
0180 preshowerSelection.add<double>("preshClusterEnergyCut", 0.0);
0181 preshowerSelection.add<double>("preshStripEnergyCut", 0.0);
0182 preshowerSelection.add<int>("preshSeededNstrip", 15);
0183 preshowerSelection.add<double>("preshCalibPlaneX", 1.0);
0184 preshowerSelection.add<double>("preshCalibPlaneY", 0.7);
0185 preshowerSelection.add<double>("preshCalibGamma", 0.024);
0186 preshowerSelection.add<double>("preshCalibMIP", 9.0E-5);
0187 preshowerSelection.add<std::string>("debugLevelES", "");
0188 desc.add<edm::ParameterSetDescription>("preshowerSelection", preshowerSelection);
0189
0190 desc.add<int>("debugLevel", 0);
0191 descriptions.add("hltEcalResonanceFilter", desc);
0192 }
0193
0194
0195 bool HLTEcalResonanceFilter::filter(edm::Event &iEvent, const edm::EventSetup &iSetup) {
0196
0197 std::unique_ptr<EBRecHitCollection> selEBRecHitCollection(new EBRecHitCollection);
0198
0199 std::unique_ptr<EERecHitCollection> selEERecHitCollection(new EERecHitCollection);
0200
0201
0202 vector<DetId> selectedEBDetIds;
0203 vector<DetId> selectedEEDetIds;
0204
0205 auto const &pTopology = iSetup.getHandle(caloTopologyRecordToken_);
0206 const CaloSubdetectorTopology *topology_eb = pTopology->getSubdetectorTopology(DetId::Ecal, EcalBarrel);
0207 const CaloSubdetectorTopology *topology_ee = pTopology->getSubdetectorTopology(DetId::Ecal, EcalEndcap);
0208
0209 auto const &geoHandle = iSetup.getHandle(caloGeometryRecordToken_);
0210 const CaloSubdetectorGeometry *geometry_es = geoHandle->getSubdetectorGeometry(DetId::Ecal, EcalPreshower);
0211 std::unique_ptr<CaloSubdetectorTopology> topology_es;
0212 if (geometry_es) {
0213 topology_es = std::make_unique<EcalPreshowerTopology>();
0214 } else {
0215 storeRecHitES_ = false;
0216 }
0217
0218
0219 edm::ESHandle<EcalChannelStatus> csHandle;
0220 if (useDBStatus_)
0221 csHandle = iSetup.getHandle(ecalChannelStatusRcdToken_);
0222 const EcalChannelStatus &channelStatus = *csHandle;
0223
0224
0225
0226 Handle<EBRecHitCollection> barrelRecHitsHandle;
0227
0228 iEvent.getByToken(barrelHitsToken_, barrelRecHitsHandle);
0229 if (!barrelRecHitsHandle.isValid()) {
0230 LogDebug("") << "AlCaEcalResonanceProducer: Error! can't get product barrel hits!" << std::endl;
0231 }
0232
0233 const EcalRecHitCollection *hitCollection_eb = barrelRecHitsHandle.product();
0234
0235 Handle<reco::BasicClusterCollection> barrelClustersHandle;
0236 iEvent.getByToken(barrelClustersToken_, barrelClustersHandle);
0237 if (!barrelClustersHandle.isValid()) {
0238 LogDebug("") << "AlCaEcalResonanceProducer: Error! can't get product barrel clusters!" << std::endl;
0239 }
0240 const reco::BasicClusterCollection *clusterCollection_eb = barrelClustersHandle.product();
0241 if (debug_ >= 1)
0242 LogDebug("") << " barrel_input_size: run " << iEvent.id().run() << " event " << iEvent.id().event()
0243 << " nhitEB: " << hitCollection_eb->size() << " nCluster: " << clusterCollection_eb->size() << endl;
0244
0245 if (doSelBarrel_) {
0246 map<int, vector<EcalRecHit> > RecHits5x5_clus;
0247 vector<int> indIsoClus;
0248 vector<int> indCandClus;
0249 vector<int> indClusSelected;
0250
0251 doSelection(EcalBarrel,
0252 clusterCollection_eb,
0253 hitCollection_eb,
0254 channelStatus,
0255 topology_eb,
0256 RecHits5x5_clus,
0257 indCandClus,
0258 indIsoClus,
0259 indClusSelected);
0260
0261
0262 for (int ind : indClusSelected) {
0263 auto it_bc3 = clusterCollection_eb->begin() + ind;
0264 const std::vector<std::pair<DetId, float> > &vid = it_bc3->hitsAndFractions();
0265 for (auto const &idItr : vid) {
0266 auto hit = hitCollection_eb->find(idItr.first);
0267 if (hit == hitCollection_eb->end())
0268 continue;
0269 selEBRecHitCollection->push_back(*hit);
0270 selectedEBDetIds.push_back(idItr.first);
0271 }
0272 }
0273
0274 if (store5x5RecHitEB_) {
0275
0276 for (int ind : indCandClus) {
0277 vector<EcalRecHit> v5x5 = RecHits5x5_clus[ind];
0278 for (auto &n : v5x5) {
0279 DetId ed = n.id();
0280 auto itdet = find(selectedEBDetIds.begin(), selectedEBDetIds.end(), ed);
0281 if (itdet == selectedEBDetIds.end()) {
0282 selectedEBDetIds.push_back(ed);
0283 selEBRecHitCollection->push_back(n);
0284 }
0285 }
0286 }
0287
0288
0289 for (int ind : indIsoClus) {
0290 auto it =
0291 find(indCandClus.begin(), indCandClus.end(), ind);
0292 if (it != indCandClus.end())
0293 continue;
0294
0295 auto it_bc3 = clusterCollection_eb->begin() + ind;
0296 DetId seedId = it_bc3->seed();
0297 std::vector<DetId> clus_v5x5 = topology_eb->getWindow(seedId, 5, 5);
0298 for (std::vector<DetId>::const_iterator idItr = clus_v5x5.begin(); idItr != clus_v5x5.end(); idItr++) {
0299 DetId ed = *idItr;
0300 auto rit = hitCollection_eb->find(ed);
0301 if (rit == hitCollection_eb->end())
0302 continue;
0303 if (!checkStatusOfEcalRecHit(channelStatus, *rit))
0304 continue;
0305 auto itdet = find(selectedEBDetIds.begin(), selectedEBDetIds.end(), ed);
0306 if (itdet == selectedEBDetIds.end()) {
0307 selectedEBDetIds.push_back(ed);
0308 selEBRecHitCollection->push_back(*rit);
0309 }
0310 }
0311 }
0312 }
0313
0314 }
0315
0316 int eb_collsize = selEBRecHitCollection->size();
0317
0318 if (debug_ >= 1)
0319 LogDebug("") << " barrel_output_size_run " << iEvent.id().run() << " event " << iEvent.id().event() << " nEBSaved "
0320 << selEBRecHitCollection->size() << std::endl;
0321
0322
0323
0324
0325 Handle<ESRecHitCollection> esRecHitsHandle;
0326 iEvent.getByToken(preshHitsToken_, esRecHitsHandle);
0327 if (!esRecHitsHandle.isValid()) {
0328 LogDebug("") << "AlCaEcalResonanceProducer: Error! can't get product esRecHit!" << std::endl;
0329 }
0330 const EcalRecHitCollection *hitCollection_es = esRecHitsHandle.product();
0331
0332 m_esrechit_map.clear();
0333 EcalRecHitCollection::const_iterator iter;
0334 for (iter = esRecHitsHandle->begin(); iter != esRecHitsHandle->end(); iter++) {
0335
0336 m_esrechit_map.insert(std::make_pair(iter->id(), *iter));
0337 }
0338
0339 m_used_strips.clear();
0340 std::unique_ptr<ESRecHitCollection> selESRecHitCollection(new ESRecHitCollection);
0341
0342 Handle<EERecHitCollection> endcapRecHitsHandle;
0343 iEvent.getByToken(endcapHitsToken_, endcapRecHitsHandle);
0344 if (!endcapRecHitsHandle.isValid()) {
0345 LogDebug("") << "AlCaEcalResonanceProducer: Error! can't get product endcap hits!" << std::endl;
0346 }
0347 const EcalRecHitCollection *hitCollection_ee = endcapRecHitsHandle.product();
0348 Handle<reco::BasicClusterCollection> endcapClustersHandle;
0349 iEvent.getByToken(endcapClustersToken_, endcapClustersHandle);
0350 if (!endcapClustersHandle.isValid()) {
0351 LogDebug("") << "AlCaEcalResonanceProducer: Error! can't get product endcap clusters!" << std::endl;
0352 }
0353 const reco::BasicClusterCollection *clusterCollection_ee = endcapClustersHandle.product();
0354 if (debug_ >= 1)
0355 LogDebug("") << " endcap_input_size: run " << iEvent.id().run() << " event " << iEvent.id().event()
0356 << " nhitEE: " << hitCollection_ee->size() << " nhitES: " << hitCollection_es->size()
0357 << " nCluster: " << clusterCollection_ee->size() << endl;
0358
0359 if (doSelEndcap_) {
0360 map<int, vector<EcalRecHit> > RecHits5x5_clus;
0361 vector<int> indIsoClus;
0362 vector<int> indCandClus;
0363 vector<int> indClusSelected;
0364
0365 doSelection(EcalEndcap,
0366 clusterCollection_ee,
0367 hitCollection_ee,
0368 channelStatus,
0369 topology_ee,
0370 RecHits5x5_clus,
0371 indCandClus,
0372 indIsoClus,
0373 indClusSelected);
0374
0375
0376 for (int ind : indClusSelected) {
0377 auto it_bc3 = clusterCollection_ee->begin() + ind;
0378 const std::vector<std::pair<DetId, float> > &vid = it_bc3->hitsAndFractions();
0379 for (auto const &idItr : vid) {
0380 auto hit = hitCollection_ee->find(idItr.first);
0381 if (hit == hitCollection_ee->end())
0382 continue;
0383 selEERecHitCollection->push_back(*hit);
0384 selectedEEDetIds.push_back(idItr.first);
0385 }
0386
0387 if (storeRecHitES_) {
0388 std::set<DetId> used_strips_before = m_used_strips;
0389 makeClusterES(it_bc3->x(), it_bc3->y(), it_bc3->z(), geometry_es, topology_es.get());
0390 auto ites = m_used_strips.begin();
0391 for (; ites != m_used_strips.end(); ++ites) {
0392 ESDetId d1 = ESDetId(*ites);
0393 auto ites2 = find(used_strips_before.begin(), used_strips_before.end(), d1);
0394 if ((ites2 == used_strips_before.end())) {
0395 auto itmap = m_esrechit_map.find(d1);
0396 selESRecHitCollection->push_back(itmap->second);
0397 }
0398 }
0399 }
0400 }
0401
0402 if (store5x5RecHitEE_) {
0403
0404 for (int ind : indCandClus) {
0405 vector<EcalRecHit> v5x5 = RecHits5x5_clus[ind];
0406 for (auto &n : v5x5) {
0407 DetId ed = n.id();
0408 auto itdet = find(selectedEEDetIds.begin(), selectedEEDetIds.end(), ed);
0409 if (itdet == selectedEEDetIds.end()) {
0410 selectedEEDetIds.push_back(ed);
0411 selEERecHitCollection->push_back(n);
0412 }
0413 }
0414 }
0415
0416
0417 for (int ind : indIsoClus) {
0418 auto it =
0419 find(indCandClus.begin(), indCandClus.end(), ind);
0420 if (it != indCandClus.end())
0421 continue;
0422
0423 auto it_bc3 = clusterCollection_ee->begin() + ind;
0424 DetId seedId = it_bc3->seed();
0425 std::vector<DetId> clus_v5x5 = topology_ee->getWindow(seedId, 5, 5);
0426 for (std::vector<DetId>::const_iterator idItr = clus_v5x5.begin(); idItr != clus_v5x5.end(); idItr++) {
0427 DetId ed = *idItr;
0428 auto rit = hitCollection_ee->find(ed);
0429 if (rit == hitCollection_ee->end())
0430 continue;
0431 if (!checkStatusOfEcalRecHit(channelStatus, *rit))
0432 continue;
0433 auto itdet = find(selectedEEDetIds.begin(), selectedEEDetIds.end(), ed);
0434 if (itdet == selectedEEDetIds.end()) {
0435 selectedEEDetIds.push_back(ed);
0436 selEERecHitCollection->push_back(*rit);
0437 }
0438 }
0439 }
0440 }
0441
0442 }
0443
0444
0445
0446 if (debug_ >= 1)
0447 LogDebug("") << " endcap_output_size run " << iEvent.id().run() << "_" << iEvent.id().event() << " nEESaved "
0448 << selEERecHitCollection->size() << " nESSaved: " << selESRecHitCollection->size() << endl;
0449
0450
0451 int ee_collsize = selEERecHitCollection->size();
0452
0453 if (eb_collsize < 2 && ee_collsize < 2)
0454 return false;
0455
0456
0457 if (doSelBarrel_) {
0458 iEvent.put(std::move(selEBRecHitCollection), BarrelHits_);
0459 }
0460 if (doSelEndcap_) {
0461 iEvent.put(std::move(selEERecHitCollection), EndcapHits_);
0462 if (storeRecHitES_) {
0463 iEvent.put(std::move(selESRecHitCollection), ESHits_);
0464 }
0465 }
0466
0467 return true;
0468 }
0469
0470 void HLTEcalResonanceFilter::doSelection(
0471 int detector,
0472 const reco::BasicClusterCollection *clusterCollection,
0473 const EcalRecHitCollection *hitCollection,
0474 const EcalChannelStatus &channelStatus,
0475 const CaloSubdetectorTopology *topology_p,
0476 map<int, vector<EcalRecHit> > &RecHits5x5_clus,
0477 vector<int> &indCandClus,
0478 vector<int> &indIsoClus,
0479 vector<int> &indClusSelected
0480 ) {
0481 vector<int> indClusPi0Candidates;
0482 if (detector == EcalBarrel && removePi0CandidatesForEta_) {
0483 for (auto it_bc = clusterCollection->begin(); it_bc != clusterCollection->end(); it_bc++) {
0484 for (auto it_bc2 = it_bc + 1; it_bc2 != clusterCollection->end(); it_bc2++) {
0485 float m_pair, pt_pair, eta_pair, phi_pair;
0486 calcPaircluster(*it_bc, *it_bc2, m_pair, pt_pair, eta_pair, phi_pair);
0487 if (m_pair > massLowPi0Cand_ && m_pair < massHighPi0Cand_) {
0488 int indtmp[2] = {int(it_bc - clusterCollection->begin()), int(it_bc2 - clusterCollection->begin())};
0489 for (int &k : indtmp) {
0490 auto it = find(indClusPi0Candidates.begin(), indClusPi0Candidates.end(), k);
0491 if (it == indClusPi0Candidates.end())
0492 indClusPi0Candidates.push_back(k);
0493 }
0494 }
0495 }
0496 }
0497 }
0498
0499 double m_minCut = seleMinvMinBarrel_;
0500 double m_maxCut = seleMinvMaxBarrel_;
0501 double ptminforIso = ptMinForIsolation_;
0502 double isoCut = seleIso_;
0503 double isoBeltdrCut = seleBeltDR_;
0504 double isoBeltdetaCut = seleBeltDeta_;
0505 double s4s9Cut = seleS4S9Gamma_;
0506 double s9s25Cut = seleS9S25Gamma_;
0507 bool store5x5 = store5x5RecHitEB_;
0508
0509 if (detector == EcalEndcap) {
0510 m_minCut = seleMinvMinEndCap_;
0511 m_maxCut = seleMinvMaxEndCap_;
0512 ptminforIso = ptMinForIsolationEndCap_;
0513 isoCut = seleIsoEndCap_;
0514 isoBeltdrCut = seleBeltDREndCap_;
0515 isoBeltdetaCut = seleBeltDetaEndCap_;
0516 s4s9Cut = seleS4S9GammaEndCap_;
0517 s9s25Cut = seleS9S25GammaEndCap_;
0518 store5x5 = store5x5RecHitEE_;
0519 }
0520
0521 map<int, bool>
0522 passShowerShape_clus;
0523
0524 for (auto it_bc = clusterCollection->begin(); it_bc != clusterCollection->end(); it_bc++) {
0525 if (detector == EcalBarrel && removePi0CandidatesForEta_) {
0526 auto it = find(indClusPi0Candidates.begin(), indClusPi0Candidates.end(), int(it_bc - clusterCollection->begin()));
0527 if (it != indClusPi0Candidates.end())
0528 continue;
0529 }
0530
0531 float en1 = it_bc->energy();
0532 float theta1 = 2. * atan(exp(-it_bc->position().eta()));
0533 float pt1 = en1 * sin(theta1);
0534
0535 int ind1 = int(it_bc - clusterCollection->begin());
0536 auto itmap = passShowerShape_clus.find(ind1);
0537 if (itmap != passShowerShape_clus.end()) {
0538 if (itmap->second == false) {
0539 continue;
0540 }
0541 }
0542
0543 for (auto it_bc2 = it_bc + 1; it_bc2 != clusterCollection->end(); it_bc2++) {
0544 if (detector == EcalBarrel && removePi0CandidatesForEta_) {
0545 auto it =
0546 find(indClusPi0Candidates.begin(), indClusPi0Candidates.end(), int(it_bc2 - clusterCollection->begin()));
0547 if (it != indClusPi0Candidates.end())
0548 continue;
0549 }
0550 float en2 = it_bc2->energy();
0551 float theta2 = 2. * atan(exp(-it_bc2->position().eta()));
0552 float pt2 = en2 * sin(theta2);
0553
0554 int ind2 = int(it_bc2 - clusterCollection->begin());
0555 auto itmap = passShowerShape_clus.find(ind2);
0556 if (itmap != passShowerShape_clus.end()) {
0557 if (itmap->second == false) {
0558 continue;
0559 }
0560 }
0561
0562 float m_pair, pt_pair, eta_pair, phi_pair;
0563 calcPaircluster(*it_bc, *it_bc2, m_pair, pt_pair, eta_pair, phi_pair);
0564
0565 float ptmin = pt1 < pt2 ? pt1 : pt2;
0566 if (detector == EcalBarrel) {
0567 if (ptmin < selePtGamma_)
0568 continue;
0569 if (pt_pair < selePtPair_)
0570 continue;
0571 } else {
0572 float etapair = fabs(eta_pair);
0573 if (etapair <= region1_EndCap_) {
0574 if (ptmin < selePtGammaEndCap_region1_ || pt_pair < selePtPairEndCap_region1_)
0575 continue;
0576 } else if (etapair <= region2_EndCap_) {
0577 if (ptmin < selePtGammaEndCap_region2_ || pt_pair < selePtPairEndCap_region2_)
0578 continue;
0579 } else {
0580 if (ptmin < selePtGammaEndCap_region3_ || pt_pair < selePtPairEndCap_region3_)
0581 continue;
0582 if (pt_pair > selePtPairMaxEndCap_region3_)
0583 continue;
0584 }
0585 }
0586
0587
0588 if (m_pair < m_minCut || m_pair > m_maxCut)
0589 continue;
0590
0591
0592 vector<int> IsoClus;
0593 IsoClus.push_back(ind1);
0594 IsoClus.push_back(ind2);
0595
0596 float Iso = 0;
0597 for (auto it_bc3 = clusterCollection->begin(); it_bc3 != clusterCollection->end(); it_bc3++) {
0598 if (it_bc3->seed() == it_bc->seed() || it_bc3->seed() == it_bc2->seed())
0599 continue;
0600 float drcl = GetDeltaR(eta_pair, it_bc3->eta(), phi_pair, it_bc3->phi());
0601 float dretacl = fabs(eta_pair - it_bc3->eta());
0602 if (drcl > isoBeltdrCut || dretacl > isoBeltdetaCut)
0603 continue;
0604 float pt3 = it_bc3->energy() * sin(it_bc3->position().theta());
0605 if (pt3 < ptminforIso)
0606 continue;
0607 Iso += pt3;
0608 int ind3 = int(it_bc3 - clusterCollection->begin());
0609 IsoClus.push_back(ind3);
0610 }
0611
0612 if (Iso / pt_pair > isoCut)
0613 continue;
0614
0615 bool failShowerShape = false;
0616 for (int n = 0; n < 2; n++) {
0617 int ind = IsoClus[n];
0618 auto it_bc3 = clusterCollection->begin() + ind;
0619 auto itmap = passShowerShape_clus.find(ind);
0620 if (itmap != passShowerShape_clus.end()) {
0621 if (itmap->second == false) {
0622 failShowerShape = true;
0623 n = 2;
0624 }
0625 } else {
0626 vector<EcalRecHit> RecHitsCluster_5x5;
0627 float res[3];
0628 calcShowerShape(*it_bc3, channelStatus, hitCollection, topology_p, store5x5, RecHitsCluster_5x5, res);
0629 float s4s9 = res[1] > 0 ? res[0] / res[1] : -1;
0630 float s9s25 = res[2] > 0 ? res[1] / res[2] : -1;
0631 bool passed = s4s9 > s4s9Cut && s9s25 > s9s25Cut;
0632 passShowerShape_clus.insert(pair<int, bool>(ind, passed));
0633 if (!passed) {
0634 failShowerShape = true;
0635 n = 2;
0636 } else {
0637 RecHits5x5_clus.insert(pair<int, vector<EcalRecHit> >(ind, RecHitsCluster_5x5));
0638 }
0639 }
0640 }
0641
0642 if (failShowerShape == true)
0643 continue;
0644
0645
0646 for (unsigned int iii = 0; iii < IsoClus.size(); iii++) {
0647 int ind = IsoClus[iii];
0648
0649 if (iii < 2) {
0650 auto it = find(indCandClus.begin(), indCandClus.end(), ind);
0651 if (it == indCandClus.end())
0652 indCandClus.push_back(ind);
0653 else
0654 continue;
0655 } else {
0656 auto it = find(indIsoClus.begin(), indIsoClus.end(), ind);
0657 if (it == indIsoClus.end())
0658 indIsoClus.push_back(ind);
0659 else
0660 continue;
0661 }
0662
0663 auto it = find(indClusSelected.begin(), indClusSelected.end(), ind);
0664 if (it != indClusSelected.end())
0665 continue;
0666 indClusSelected.push_back(ind);
0667 }
0668 }
0669 }
0670 }
0671
0672 void HLTEcalResonanceFilter::convxtalid(Int_t &nphi, Int_t &neta) {
0673
0674
0675
0676
0677 if (neta > 0)
0678 neta -= 1;
0679 if (nphi > 359)
0680 nphi = nphi - 360;
0681
0682
0683 if (nphi > 359 || nphi < 0 || neta < -85 || neta > 84) {
0684 LogError("") << " unexpected fatal error in HLTEcalResonanceFilter::convxtalid " << nphi << " " << neta << " "
0685 << std::endl;
0686 }
0687 }
0688
0689 int HLTEcalResonanceFilter::diff_neta_s(Int_t neta1, Int_t neta2) {
0690 Int_t mdiff;
0691 mdiff = (neta1 - neta2);
0692 return mdiff;
0693 }
0694
0695
0696 int HLTEcalResonanceFilter::diff_nphi_s(Int_t nphi1, Int_t nphi2) {
0697 Int_t mdiff;
0698 if (std::abs(nphi1 - nphi2) < (360 - std::abs(nphi1 - nphi2))) {
0699 mdiff = nphi1 - nphi2;
0700 } else {
0701 mdiff = 360 - std::abs(nphi1 - nphi2);
0702 if (nphi1 > nphi2)
0703 mdiff = -mdiff;
0704 }
0705 return mdiff;
0706 }
0707
0708 void HLTEcalResonanceFilter::calcShowerShape(const reco::BasicCluster &bc,
0709 const EcalChannelStatus &channelStatus,
0710 const EcalRecHitCollection *recHits,
0711 const CaloSubdetectorTopology *topology_p,
0712 bool calc5x5,
0713 vector<EcalRecHit> &rechit5x5,
0714 float res[]) {
0715 const std::vector<std::pair<DetId, float> > &vid = bc.hitsAndFractions();
0716
0717 float e2x2[4] = {0};
0718 float e3x3 = 0;
0719 float e5x5 = 0;
0720 int seedx;
0721 int seedy;
0722 DetId seedId = bc.seed();
0723
0724 bool InBarrel = true;
0725 if (seedId.subdetId() == EcalBarrel) {
0726 EBDetId ebd = EBDetId(seedId);
0727 seedx = ebd.ieta();
0728 seedy = ebd.iphi();
0729 convxtalid(seedy, seedx);
0730 } else {
0731 InBarrel = false;
0732 EEDetId eed = EEDetId(seedId);
0733 seedx = eed.ix();
0734 seedy = eed.iy();
0735 }
0736 int x, y, dx, dy;
0737
0738 if (calc5x5) {
0739 rechit5x5.clear();
0740 std::vector<DetId> clus_v5x5;
0741 clus_v5x5 = topology_p->getWindow(seedId, 5, 5);
0742
0743 for (std::vector<DetId>::const_iterator idItr = clus_v5x5.begin(); idItr != clus_v5x5.end(); idItr++) {
0744 DetId ed = *idItr;
0745 if (InBarrel == true) {
0746 EBDetId ebd = EBDetId(ed);
0747 x = ebd.ieta();
0748 y = ebd.iphi();
0749 convxtalid(y, x);
0750 dx = diff_neta_s(x, seedx);
0751 dy = diff_nphi_s(y, seedy);
0752 } else {
0753 EEDetId eed = EEDetId(ed);
0754 x = eed.ix();
0755 y = eed.iy();
0756 dx = x - seedx;
0757 dy = y - seedy;
0758 }
0759 auto rit = recHits->find(ed);
0760 if (rit == recHits->end())
0761 continue;
0762 if (!checkStatusOfEcalRecHit(channelStatus, *rit))
0763 continue;
0764
0765 float energy = (*rit).energy();
0766 e5x5 += energy;
0767
0768 auto idItrF = std::find(vid.begin(), vid.end(), std::make_pair(ed, 1.0f));
0769 if (idItrF == vid.end()) {
0770 rechit5x5.push_back(*rit);
0771 } else {
0772 if (std::abs(dx) <= 1 && std::abs(dy) <= 1) {
0773 if (dx <= 0 && dy <= 0)
0774 e2x2[0] += energy;
0775 if (dx >= 0 && dy <= 0)
0776 e2x2[1] += energy;
0777 if (dx <= 0 && dy >= 0)
0778 e2x2[2] += energy;
0779 if (dx >= 0 && dy >= 0)
0780 e2x2[3] += energy;
0781 e3x3 += energy;
0782 }
0783 }
0784 }
0785
0786 } else {
0787 for (auto const &idItr : vid) {
0788 DetId ed = idItr.first;
0789 if (InBarrel == true) {
0790 EBDetId ebd = EBDetId(ed);
0791 x = ebd.ieta();
0792 y = ebd.iphi();
0793 convxtalid(y, x);
0794 dx = diff_neta_s(x, seedx);
0795 dy = diff_nphi_s(y, seedy);
0796 } else {
0797 EEDetId eed = EEDetId(ed);
0798 x = eed.ix();
0799 y = eed.iy();
0800 dx = x - seedx;
0801 dy = y - seedy;
0802 }
0803 auto rit = recHits->find(ed);
0804 if (rit == recHits->end()) {
0805 continue;
0806 }
0807
0808 float energy = (*rit).energy();
0809 if (std::abs(dx) <= 1 && std::abs(dy) <= 1) {
0810 if (dx <= 0 && dy <= 0)
0811 e2x2[0] += energy;
0812 if (dx >= 0 && dy <= 0)
0813 e2x2[1] += energy;
0814 if (dx <= 0 && dy >= 0)
0815 e2x2[2] += energy;
0816 if (dx >= 0 && dy >= 0)
0817 e2x2[3] += energy;
0818 e3x3 += energy;
0819 }
0820 }
0821 e5x5 = e3x3;
0822 }
0823
0824
0825 res[0] = *max_element(e2x2, e2x2 + 4);
0826 res[1] = e3x3;
0827 res[2] = e5x5;
0828 }
0829
0830 void HLTEcalResonanceFilter::calcPaircluster(const reco::BasicCluster &bc1,
0831 const reco::BasicCluster &bc2,
0832 float &m_pair,
0833 float &pt_pair,
0834 float &eta_pair,
0835 float &phi_pair) {
0836 float theta1 = 2. * atan(exp(-bc1.eta()));
0837 float en1 = bc1.energy();
0838 float pt1 = en1 * sin(theta1);
0839 TLorentzVector v1(pt1 * cos(bc1.phi()), pt1 * sin(bc1.phi()), en1 * cos(theta1), en1);
0840
0841 float theta2 = 2. * atan(exp(-bc2.eta()));
0842 float en2 = bc2.energy();
0843 float pt2 = en2 * sin(theta2);
0844 TLorentzVector v2(pt2 * cos(bc2.phi()), pt2 * sin(bc2.phi()), en2 * cos(theta2), en2);
0845
0846 TLorentzVector v = v1 + v2;
0847
0848 m_pair = v.M();
0849 pt_pair = v.Pt();
0850 eta_pair = v.Eta();
0851 phi_pair = v.Phi();
0852 }
0853
0854 bool HLTEcalResonanceFilter::checkStatusOfEcalRecHit(const EcalChannelStatus &channelStatus, const EcalRecHit &rh) {
0855 if (useRecoFlag_) {
0856 int flag = rh.recoFlag();
0857 if (flagLevelRecHitsToUse_ == 0) {
0858 if (flag != 0)
0859 return false;
0860 } else if (flagLevelRecHitsToUse_ == 1) {
0861 if (flag != 0 && flag != 4)
0862 return false;
0863 } else if (flagLevelRecHitsToUse_ == 2) {
0864 if (flag != 0 && flag != 4 && flag != 6 && flag != 7)
0865 return false;
0866 }
0867 }
0868 if (useDBStatus_) {
0869 int status = int(channelStatus[rh.id().rawId()].getStatusCode());
0870 if (status > statusLevelRecHitsToUse_)
0871 return false;
0872 }
0873
0874 return true;
0875 }
0876
0877 float HLTEcalResonanceFilter::DeltaPhi(float phi1, float phi2) {
0878 float diff = fabs(phi2 - phi1);
0879
0880 while (diff > acos(-1))
0881 diff -= 2 * acos(-1);
0882 while (diff <= -acos(-1))
0883 diff += 2 * acos(-1);
0884
0885 return diff;
0886 }
0887
0888 float HLTEcalResonanceFilter::GetDeltaR(float eta1, float eta2, float phi1, float phi2) {
0889 return sqrt((eta1 - eta2) * (eta1 - eta2) + DeltaPhi(phi1, phi2) * DeltaPhi(phi1, phi2));
0890 }
0891
0892 void HLTEcalResonanceFilter::makeClusterES(
0893 float x, float y, float z, const CaloSubdetectorGeometry *geometry_es, const CaloSubdetectorTopology *topology_es) {
0894
0895 const GlobalPoint point(x, y, z);
0896 DetId tmp1 = (dynamic_cast<const EcalPreshowerGeometry *>(geometry_es))->getClosestCellInPlane(point, 1);
0897 DetId tmp2 = (dynamic_cast<const EcalPreshowerGeometry *>(geometry_es))->getClosestCellInPlane(point, 2);
0898 ESDetId strip1 = (tmp1 == DetId(0)) ? ESDetId(0) : ESDetId(tmp1);
0899 ESDetId strip2 = (tmp2 == DetId(0)) ? ESDetId(0) : ESDetId(tmp2);
0900
0901
0902 for (int i2 = 0; i2 < preshNclust_; i2++) {
0903 reco::PreshowerCluster cl1 =
0904 presh_algo_->makeOneCluster(strip1, &m_used_strips, &m_esrechit_map, geometry_es, topology_es);
0905 reco::PreshowerCluster cl2 =
0906 presh_algo_->makeOneCluster(strip2, &m_used_strips, &m_esrechit_map, geometry_es, topology_es);
0907 }
0908 }
0909
0910
0911 #include "FWCore/Framework/interface/MakerMacros.h"
0912 DEFINE_FWK_MODULE(HLTEcalResonanceFilter);