Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:18:40 

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