Back to home page

Project CMSSW displayed by LXR
 
 

    

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     ///for barrel selection
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     ///for endcap selection
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   ///for storing rechits ES for each selected EE clusters.
0093   storeRecHitES_ = iConfig.getParameter<bool>("storeRecHitES");
0094   if (storeRecHitES_) {
0095     edm::ParameterSet preshowerSelection = iConfig.getParameter<edm::ParameterSet>("preshowerSelection");
0096 
0097     // maximum number of matched ES clusters (in each ES layer) to each BC
0098     preshNclust_ = preshowerSelection.getParameter<int>("preshNclust");
0099     // min energy of ES clusters
0100     preshClustECut = preshowerSelection.getParameter<double>("preshClusterEnergyCut");
0101     // algo params
0102     float preshStripECut = preshowerSelection.getParameter<double>("preshStripEnergyCut");
0103     int preshSeededNst = preshowerSelection.getParameter<int>("preshSeededNstrip");
0104     // calibration parameters:
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     // ES algo constructor:
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", "");  // *** This is not needed and shoul be better removed !
0188   desc.add<edm::ParameterSetDescription>("preshowerSelection", preshowerSelection);
0189 
0190   desc.add<int>("debugLevel", 0);
0191   descriptions.add("hltEcalResonanceFilter", desc);
0192 }
0193 
0194 // ------------ method called to produce the data  ------------
0195 bool HLTEcalResonanceFilter::filter(edm::Event &iEvent, const edm::EventSetup &iSetup) {
0196   //Create empty output collections
0197   std::unique_ptr<EBRecHitCollection> selEBRecHitCollection(new EBRecHitCollection);
0198   //Create empty output collections
0199   std::unique_ptr<EERecHitCollection> selEERecHitCollection(new EERecHitCollection);
0200 
0201   ////all selected..
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;  ///if no preshower
0216   }
0217 
0218   ///get status from DB
0219   edm::ESHandle<EcalChannelStatus> csHandle;
0220   if (useDBStatus_)
0221     csHandle = iSetup.getHandle(ecalChannelStatusRcdToken_);
0222   const EcalChannelStatus &channelStatus = *csHandle;
0223 
0224   ///==============Start to process barrel part==================///
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;  ///5x5 for selected pairs
0247     vector<int> indIsoClus;                         /// Iso cluster all , 5x5 rechit not yet done
0248     vector<int> indCandClus;                        ///good cluster all ,  5x5 rechit done already during the loop
0249     vector<int> indClusSelected;                    /// saved so far, all
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     ///Now save all rechits in the selected clusters
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;  //this should not happen.
0269         selEBRecHitCollection->push_back(*hit);
0270         selectedEBDetIds.push_back(idItr.first);
0271       }
0272     }
0273 
0274     if (store5x5RecHitEB_) {
0275       ///stroe 5x5 of good clusters, 5x5 arleady got
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       ///store 5x5 of Iso clusters, need to getWindow of 5x5
0289       for (int ind : indIsoClus) {
0290         auto it =
0291             find(indCandClus.begin(), indCandClus.end(), ind);  ///check if already saved in the good cluster vector
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   }  // end of selection for eta/pi0->gg in the barrel
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   ///==============End of  barrel ==================///
0322 
0323   //===============Start of Endcap =================/////
0324   ///get preshower rechits
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   // make a map of rechits:
0332   m_esrechit_map.clear();
0333   EcalRecHitCollection::const_iterator iter;
0334   for (iter = esRecHitsHandle->begin(); iter != esRecHitsHandle->end(); iter++) {
0335     //Make the map of DetID, EcalRecHit pairs
0336     m_esrechit_map.insert(std::make_pair(iter->id(), *iter));
0337   }
0338   // The set of used DetID's for a given event:
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;  ///5x5 for selected pairs
0361     vector<int> indIsoClus;                         /// Iso cluster all , 5x5 rechit not yet done
0362     vector<int> indCandClus;                        ///good cluster all ,  5x5 rechit done already during the loop
0363     vector<int> indClusSelected;                    /// saved so far, all
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     ///Now save all rechits in the selected clusters
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;  //this should not happen.
0383         selEERecHitCollection->push_back(*hit);
0384         selectedEEDetIds.push_back(idItr.first);
0385       }
0386       ///save preshower rechits
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       ///store 5x5 of good clusters, 5x5 arleady got
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       ///store 5x5 of Iso clusters, need to getWindow of 5x5
0417       for (int ind : indIsoClus) {
0418         auto it =
0419             find(indCandClus.begin(), indCandClus.end(), ind);  ///check if already saved in the good cluster vector
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   }  // end of selection for eta/pi0->gg in the endcap
0443 
0444   ////==============End of endcap ===============///
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   //Put selected information in the event
0451   int ee_collsize = selEERecHitCollection->size();
0452 
0453   if (eb_collsize < 2 && ee_collsize < 2)
0454     return false;
0455 
0456   ////Now put into events selected rechits.
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,     ///good cluster all ,  5x5 rechit done already during the loop
0478     vector<int> &indIsoClus,      /// Iso cluster all , 5x5 rechit not yet done
0479     vector<int> &indClusSelected  /// saved so far, all
0480 ) {
0481   vector<int> indClusPi0Candidates;  ///those clusters identified as pi0s
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     }  //end of loop over finding pi0's clusters
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;  //if this cluster passed showershape cut, so no need to compute the quantity again for each loop
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       /// pt & Pt pair Cut
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       /// Mass window Cut
0588       if (m_pair < m_minCut || m_pair > m_maxCut)
0589         continue;
0590 
0591       //// Loop on cluster to measure isolation:
0592       vector<int> IsoClus;
0593       IsoClus.push_back(ind1);  //first two are good clusters
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());  /// remember which Iso cluster used
0609         IsoClus.push_back(ind3);
0610       }
0611       /// Isolation cut
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;  //exit the loop
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;  //exit the loop
0636           } else {
0637             RecHits5x5_clus.insert(pair<int, vector<EcalRecHit> >(ind, RecHitsCluster_5x5));
0638           }
0639         }
0640       }
0641 
0642       if (failShowerShape == true)
0643         continue;  //if any of the two clusters fail shower shape
0644 
0645       ///Save two good clusters' index( plus Iso cluster )  if not yet saved
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);  ///good cluster
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);  //iso cluster
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   // Barrel only
0674   // Output nphi 0...359; neta 0...84; nside=+1 (for eta>0), or 0 (for eta<0).
0675   // neta will be [-85,-1] , or [0,84], the minus sign indicates the z<0 side.
0676 
0677   if (neta > 0)
0678     neta -= 1;
0679   if (nphi > 359)
0680     nphi = nphi - 360;
0681 
0682   // final check
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 }  //end of convxtalid
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 // Calculate the distance in xtals taking into account the periodicity of the Barrel
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));  ///has to add "f", make it float
0769       if (idItrF == vid.end()) {  ///only store those not belonging to this cluster
0770         rechit5x5.push_back(*rit);
0771       } else {  /// S4, S9 are defined inside the cluster, the same as below.
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;  ///if not asked to calculte 5x5, then just make e5x5 = e3x3
0822   }
0823 
0824   ///e2x2
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_) {  ///from recoFlag()
0856     int flag = rh.recoFlag();
0857     if (flagLevelRecHitsToUse_ == 0) {  ///good
0858       if (flag != 0)
0859         return false;
0860     } else if (flagLevelRecHitsToUse_ == 1) {  ///good || PoorCalib
0861       if (flag != 0 && flag != 4)
0862         return false;
0863     } else if (flagLevelRecHitsToUse_ == 2) {  ///good || PoorCalib || LeadingEdgeRecovered || kNeighboursRecovered,
0864       if (flag != 0 && flag != 4 && flag != 6 && flag != 7)
0865         return false;
0866     }
0867   }
0868   if (useDBStatus_) {  //// from DB
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   ///get assosicated ES clusters of this endcap cluster
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   // Get ES clusters (found by the PreshSeeded algorithm) associated with a given EE cluster.
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   }  // end of cycle over ES clusters
0908 }
0909 
0910 // declare this class as a framework plugin
0911 #include "FWCore/Framework/interface/MakerMacros.h"
0912 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!