Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-11-06 06:06:43 

0001 #include "HLTRegionalEcalResonanceFilter.h"
0002 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0003 #include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
0004 
0005 #include "DataFormats/Math/interface/Vector3D.h"  // to use math::XYZVector
0006 #include "DataFormats/Math/interface/deltaPhi.h"  // use already existing deltaPhi function to return deltaPhi in [-pi, pi]
0007 
0008 using namespace std;
0009 using namespace edm;
0010 
0011 HLTRegionalEcalResonanceFilter::HLTRegionalEcalResonanceFilter(const edm::ParameterSet &iConfig)
0012     : caloTopologyRecordToken_(esConsumes()),
0013       ecalChannelStatusRcdToken_(esConsumes()),
0014       caloGeometryRecordToken_(esConsumes()) {
0015   barrelHits_ = iConfig.getParameter<edm::InputTag>("barrelHits");
0016   barrelClusters_ = iConfig.getParameter<edm::InputTag>("barrelClusters");
0017   barrelHitsToken_ = consumes<EBRecHitCollection>(barrelHits_);
0018   barrelClustersToken_ = consumes<reco::BasicClusterCollection>(barrelClusters_);
0019 
0020   endcapHits_ = iConfig.getParameter<edm::InputTag>("endcapHits");
0021   endcapClusters_ = iConfig.getParameter<edm::InputTag>("endcapClusters");
0022   endcapHitsToken_ = consumes<EERecHitCollection>(endcapHits_);
0023   endcapClustersToken_ = consumes<reco::BasicClusterCollection>(endcapClusters_);
0024 
0025   store5x5RecHitEB_ = false;
0026 
0027   doSelBarrel_ = iConfig.getParameter<bool>("doSelBarrel");
0028 
0029   if (doSelBarrel_) {
0030     edm::ParameterSet barrelSelection = iConfig.getParameter<edm::ParameterSet>("barrelSelection");
0031 
0032     ///---------------------------BARREL SELECTION-----------------------------------
0033     // EB region 1
0034     region1_Barrel_ =
0035         barrelSelection.getParameter<double>("region1_Barrel");  //eta dividing between region 1 and region 2
0036     selePtGammaBarrel_region1_ = barrelSelection.getParameter<double>("selePtGammaBarrel_region1");
0037     selePtPairBarrel_region1_ = barrelSelection.getParameter<double>("selePtPairBarrel_region1");
0038     seleS4S9GammaBarrel_region1_ = barrelSelection.getParameter<double>("seleS4S9GammaBarrel_region1");
0039     seleIsoBarrel_region1_ = barrelSelection.getParameter<double>("seleIsoBarrel_region1");
0040     seleNxtalBarrel_region1_ = barrelSelection.getParameter<uint32_t>("seleNxtalBarrel_region1");
0041 
0042     // EB region 2
0043     selePtGammaBarrel_region2_ = barrelSelection.getParameter<double>("selePtGammaBarrel_region2");
0044     selePtPairBarrel_region2_ = barrelSelection.getParameter<double>("selePtPairBarrel_region2");
0045     seleS4S9GammaBarrel_region2_ = barrelSelection.getParameter<double>("seleS4S9GammaBarrel_region2");
0046     seleIsoBarrel_region2_ = barrelSelection.getParameter<double>("seleIsoBarrel_region2");
0047     seleNxtalBarrel_region2_ = barrelSelection.getParameter<uint32_t>("seleNxtalBarrel_region2");
0048 
0049     // other
0050     seleS9S25Gamma_ = barrelSelection.getParameter<double>("seleS9S25Gamma");
0051 
0052     //mass window
0053     seleMinvMaxBarrel_ = barrelSelection.getParameter<double>("seleMinvMaxBarrel");
0054     seleMinvMinBarrel_ = barrelSelection.getParameter<double>("seleMinvMinBarrel");
0055 
0056     // remove pi0 candidates for eta dataset
0057     removePi0CandidatesForEta_ = barrelSelection.getParameter<bool>("removePi0CandidatesForEta");
0058     if (removePi0CandidatesForEta_) {
0059       massLowPi0Cand_ = barrelSelection.getParameter<double>("massLowPi0Cand");
0060       massHighPi0Cand_ = barrelSelection.getParameter<double>("massHighPi0Cand");
0061     }
0062 
0063     // EB Isolation configuration
0064     ptMinForIsolation_ = barrelSelection.getParameter<double>("ptMinForIsolation");
0065     seleBeltDR_ = barrelSelection.getParameter<double>("seleBeltDR");
0066     seleBeltDeta_ = barrelSelection.getParameter<double>("seleBeltDeta");
0067 
0068     // EB storage and collection
0069     store5x5RecHitEB_ = barrelSelection.getParameter<bool>("store5x5RecHitEB");
0070     BarrelHits_ = barrelSelection.getParameter<string>("barrelHitCollection");
0071 
0072     // selePtGamma_ = barrelSelection.getParameter<double> ("selePtGamma");  // old non-region filter
0073     // selePtPair_ = barrelSelection.getParameter<double> ("selePtPair");   // old non-region filter
0074     // seleS4S9Gamma_ = barrelSelection.getParameter<double> ("seleS4S9Gamma");  //old non-region filter
0075     // seleIso_ = barrelSelection.getParameter<double> ("seleIso");   // old non-region filter
0076 
0077     produces<EBRecHitCollection>(BarrelHits_);
0078   }
0079 
0080   doSelEndcap_ = iConfig.getParameter<bool>("doSelEndcap");
0081   if (doSelEndcap_) {
0082     edm::ParameterSet endcapSelection = iConfig.getParameter<edm::ParameterSet>("endcapSelection");
0083 
0084     ///---------------------------ENDCAP SELECTION-----------------------------
0085     seleMinvMaxEndCap_ = endcapSelection.getParameter<double>("seleMinvMaxEndCap");
0086     seleMinvMinEndCap_ = endcapSelection.getParameter<double>("seleMinvMinEndCap");
0087 
0088     // EE region 1
0089     region1_EndCap_ =
0090         endcapSelection.getParameter<double>("region1_EndCap");  //eta dividing between region 1 and region 2
0091     selePtGammaEndCap_region1_ = endcapSelection.getParameter<double>("selePtGammaEndCap_region1");
0092     selePtPairEndCap_region1_ = endcapSelection.getParameter<double>("selePtPairEndCap_region1");
0093     seleS4S9GammaEndCap_region1_ = endcapSelection.getParameter<double>("seleS4S9GammaEndCap_region1");
0094     seleIsoEndCap_region1_ = endcapSelection.getParameter<double>("seleIsoEndCap_region1");
0095     seleNxtalEndCap_region1_ = endcapSelection.getParameter<uint32_t>("seleNxtalEndCap_region1");
0096 
0097     // EE region 2
0098     region2_EndCap_ =
0099         endcapSelection.getParameter<double>("region2_EndCap");  //eta dividing between region 2 and region 3
0100     selePtGammaEndCap_region2_ = endcapSelection.getParameter<double>("selePtGammaEndCap_region2");
0101     selePtPairEndCap_region2_ = endcapSelection.getParameter<double>("selePtPairEndCap_region2");
0102     seleS4S9GammaEndCap_region2_ = endcapSelection.getParameter<double>("seleS4S9GammaEndCap_region2");
0103     seleIsoEndCap_region2_ = endcapSelection.getParameter<double>("seleIsoEndCap_region2");
0104     seleNxtalEndCap_region2_ = endcapSelection.getParameter<uint32_t>("seleNxtalEndCap_region2");
0105 
0106     // EE region 3 (available but not yet used)
0107     selePtGammaEndCap_region3_ = endcapSelection.getParameter<double>("selePtGammaEndCap_region3");
0108     selePtPairEndCap_region3_ = endcapSelection.getParameter<double>("selePtPairEndCap_region3");
0109     selePtPairMaxEndCap_region3_ = endcapSelection.getParameter<double>("selePtPairMaxEndCap_region3");
0110     seleS4S9GammaEndCap_region3_ = endcapSelection.getParameter<double>("seleS4S9GammaEndCap_region3");
0111     seleIsoEndCap_region3_ = endcapSelection.getParameter<double>("seleIsoEndCap_region3");
0112     seleNxtalEndCap_region3_ = endcapSelection.getParameter<uint32_t>("seleNxtalEndCap_region3");
0113 
0114     seleS9S25GammaEndCap_ = endcapSelection.getParameter<double>("seleS9S25GammaEndCap");
0115 
0116     // isolation belt and size configuration
0117     ptMinForIsolationEndCap_ = endcapSelection.getParameter<double>("ptMinForIsolationEndCap");
0118     seleBeltDREndCap_ = endcapSelection.getParameter<double>("seleBeltDREndCap");
0119     seleBeltDetaEndCap_ = endcapSelection.getParameter<double>("seleBeltDetaEndCap");
0120 
0121     // EE storage and collections
0122     store5x5RecHitEE_ = endcapSelection.getParameter<bool>("store5x5RecHitEE");
0123     EndcapHits_ = endcapSelection.getParameter<string>("endcapHitCollection");
0124 
0125     //seleS4S9GammaEndCap_ = endcapSelection.getParameter<double> ("seleS4S9GammaEndCap");   // old non-region filter
0126     //seleIsoEndCap_ = endcapSelection.getParameter<double> ("seleIsoEndCap");   // old non-region filter
0127 
0128     produces<EERecHitCollection>(EndcapHits_);
0129   }
0130 
0131   useRecoFlag_ = iConfig.getParameter<bool>("useRecoFlag");
0132   flagLevelRecHitsToUse_ = iConfig.getParameter<int>("flagLevelRecHitsToUse");
0133 
0134   useDBStatus_ = iConfig.getParameter<bool>("useDBStatus");
0135   statusLevelRecHitsToUse_ = iConfig.getParameter<int>("statusLevelRecHitsToUse");
0136 
0137   preshHitProducer_ = iConfig.getParameter<edm::InputTag>("preshRecHitProducer");
0138   preshHitsToken_ = consumes<EBRecHitCollection>(preshHitProducer_);
0139 
0140   ///for storing rechits ES for each selected EE clusters.
0141   storeRecHitES_ = iConfig.getParameter<bool>("storeRecHitES");
0142   if (storeRecHitES_) {
0143     edm::ParameterSet preshowerSelection = iConfig.getParameter<edm::ParameterSet>("preshowerSelection");
0144 
0145     // maximum number of matched ES clusters (in each ES layer) to each BC
0146     preshNclust_ = preshowerSelection.getParameter<int>("preshNclust");
0147     // min energy of ES clusters
0148     preshClustECut = preshowerSelection.getParameter<double>("preshClusterEnergyCut");
0149     // algo params
0150     float preshStripECut = preshowerSelection.getParameter<double>("preshStripEnergyCut");
0151     int preshSeededNst = preshowerSelection.getParameter<int>("preshSeededNstrip");
0152     // calibration parameters:
0153     calib_planeX_ = preshowerSelection.getParameter<double>("preshCalibPlaneX");
0154     calib_planeY_ = preshowerSelection.getParameter<double>("preshCalibPlaneY");
0155     gamma_ = preshowerSelection.getParameter<double>("preshCalibGamma");
0156     mip_ = preshowerSelection.getParameter<double>("preshCalibMIP");
0157 
0158     // ES algo constructor:
0159     presh_algo_ = new PreshowerClusterAlgo(preshStripECut, preshClustECut, preshSeededNst);
0160 
0161     ESHits_ = preshowerSelection.getParameter<std::string>("ESCollection");
0162     produces<ESRecHitCollection>(ESHits_);
0163   }
0164 
0165   debug_ = iConfig.getParameter<int>("debugLevel");
0166 }
0167 
0168 HLTRegionalEcalResonanceFilter::~HLTRegionalEcalResonanceFilter() {
0169   if (storeRecHitES_) {
0170     delete presh_algo_;
0171   }
0172 }
0173 
0174 void HLTRegionalEcalResonanceFilter::fillDescriptions(edm::ConfigurationDescriptions &descriptions) {
0175   edm::ParameterSetDescription desc;
0176   desc.add<edm::InputTag>("barrelHits", edm::InputTag("hltEcalRegionalPi0EtaRecHit", "EcalRecHitsEB"));
0177   desc.add<edm::InputTag>("endcapHits", edm::InputTag("hltEcalRegionalPi0EtaRecHit", "EcalRecHitsEE"));
0178   desc.add<edm::InputTag>("preshRecHitProducer", edm::InputTag("hltEcalRegionalPi0EtaRecHit", "EcalRecHitsES"));
0179   desc.add<edm::InputTag>("barrelClusters", edm::InputTag("hltSimple3x3Clusters", "Simple3x3ClustersBarrel"));
0180   desc.add<edm::InputTag>("endcapClusters", edm::InputTag("hltSimple3x3Clusters", "Simple3x3ClustersEndcap"));
0181   desc.add<bool>("useRecoFlag", false);
0182   desc.add<int>("flagLevelRecHitsToUse", 1);
0183   desc.add<bool>("useDBStatus", true);
0184   desc.add<int>("statusLevelRecHitsToUse", 1);
0185 
0186   //----------------------BARREL CONFIGURATION-----------------
0187   desc.add<bool>("doSelBarrel", true);
0188   edm::ParameterSetDescription barrelSelection;
0189 
0190   //EB region 1
0191   barrelSelection.add<double>("region1_Barrel", 1.0);  //separator between barrel region 1 and region 2
0192   barrelSelection.add<double>("selePtGammaBarrel_region1", 1.0);
0193   barrelSelection.add<double>("selePtPairBarrel_region1", 2.0);
0194   barrelSelection.add<double>("seleIsoBarrel_region1", 0.5);
0195   barrelSelection.add<uint32_t>("seleNxtalBarrel_region1", 6);
0196   barrelSelection.add<double>("seleS4S9GammaBarrel_region1", 0.83);
0197 
0198   //EB region 2
0199   barrelSelection.add<double>("selePtGammaBarrel_region2", 1.0);
0200   barrelSelection.add<double>("selePtPairBarrel_region2", 2.0);
0201   barrelSelection.add<double>("seleIsoBarrel_region2", 0.5);
0202   barrelSelection.add<uint32_t>("seleNxtalBarrel_region2", 6);
0203   barrelSelection.add<double>("seleS4S9GammaBarrel_region2", 0.83);
0204 
0205   //EB Isolation configuration
0206   barrelSelection.add<double>("ptMinForIsolation", 1.0);
0207   barrelSelection.add<double>("seleBeltDR", 0.2);
0208   barrelSelection.add<double>("seleBeltDeta", 0.05);
0209 
0210   //other parameters
0211   barrelSelection.add<double>("seleMinvMaxBarrel", 0.22);
0212   barrelSelection.add<double>("seleMinvMinBarrel", 0.06);
0213   barrelSelection.add<bool>("removePi0CandidatesForEta", false);
0214   barrelSelection.add<double>("massLowPi0Cand", 0.104);
0215   barrelSelection.add<double>("massHighPi0Cand", 0.163);
0216   barrelSelection.add<double>("seleS9S25Gamma", 0.);
0217 
0218   //collections and storage
0219   barrelSelection.add<bool>("store5x5RecHitEB", false);
0220   barrelSelection.add<std::string>("barrelHitCollection", "pi0EcalRecHitsEB");
0221   desc.add<edm::ParameterSetDescription>("barrelSelection", barrelSelection);
0222 
0223   //barrelSelection.add<double>("selePtGamma",1.); //old non-region
0224   //barrelSelection.add<double>("selePtPair",2.); //old non-regional
0225   //barrelSelection.add<double>("seleIso",0.5); //old non-regional
0226   //barrelSelection.add<double>("seleS4S9Gamma",0.83); //old non-regional
0227 
0228   //----------------------ENDCAP CONFIGURATION-----------------
0229 
0230   desc.add<bool>("doSelEndcap", true);
0231   edm::ParameterSetDescription endcapSelection;
0232   // Mass Cuts
0233   endcapSelection.add<double>("seleMinvMaxEndCap", 0.3);
0234   endcapSelection.add<double>("seleMinvMinEndCap", 0.05);
0235 
0236   // EE region 1
0237   endcapSelection.add<double>("region1_EndCap", 2.0);  // eta division between endcap region 1 and 2
0238   endcapSelection.add<double>("selePtGammaEndCap_region1", 0.8);
0239   endcapSelection.add<double>("selePtPairEndCap_region1", 3.0);
0240   endcapSelection.add<double>("seleS4S9GammaEndCap_region1", 0.9);
0241   endcapSelection.add<double>("seleIsoEndCap_region1", 0.5);
0242   endcapSelection.add<uint32_t>("seleNxtalEndCap_region1", 6);
0243 
0244   // EE region 2
0245   endcapSelection.add<double>("region2_EndCap", 2.5);  // eta division between endcap region 2 and 3
0246   endcapSelection.add<double>("selePtGammaEndCap_region2", 0.5);
0247   endcapSelection.add<double>("selePtPairEndCap_region2", 2.0);
0248   endcapSelection.add<double>("seleS4S9GammaEndCap_region2", 0.9);
0249   endcapSelection.add<double>("seleIsoEndCap_region2", 0.5);
0250   endcapSelection.add<uint32_t>("seleNxtalEndCap_region2", 6);
0251 
0252   // EE region 3
0253   endcapSelection.add<double>("selePtGammaEndCap_region3", 0.3);
0254   endcapSelection.add<double>("selePtPairEndCap_region3", 1.2);
0255   endcapSelection.add<double>("selePtPairMaxEndCap_region3", 2.5);
0256   endcapSelection.add<double>("seleS4S9GammaEndCap_region3", 0.9);
0257   endcapSelection.add<double>("seleIsoEndCap_region3", 0.5);
0258   endcapSelection.add<uint32_t>("seleNxtalEndCap_region3", 6);
0259 
0260   // other
0261   endcapSelection.add<double>("seleS9S25GammaEndCap", 0.);
0262 
0263   // isolation configuration for endcap
0264   endcapSelection.add<double>("ptMinForIsolationEndCap", 0.5);
0265   endcapSelection.add<double>("seleBeltDREndCap", 0.2);
0266   endcapSelection.add<double>("seleBeltDetaEndCap", 0.05);
0267 
0268   // collections and storage
0269   endcapSelection.add<bool>("store5x5RecHitEE", false);
0270   endcapSelection.add<std::string>("endcapHitCollection", "pi0EcalRecHitsEE");
0271   desc.add<edm::ParameterSetDescription>("endcapSelection", endcapSelection);
0272 
0273   //endcapSelection.add<double>("seleS4S9GammaEndCap",0.9); //old non-region filter
0274 
0275   //-----------------------------------------------------------
0276 
0277   desc.add<bool>("storeRecHitES", true);
0278   edm::ParameterSetDescription preshowerSelection;
0279   preshowerSelection.add<std::string>("ESCollection", "pi0EcalRecHitsES");
0280   preshowerSelection.add<int>("preshNclust", 4);
0281   preshowerSelection.add<double>("preshClusterEnergyCut", 0.0);
0282   preshowerSelection.add<double>("preshStripEnergyCut", 0.0);
0283   preshowerSelection.add<int>("preshSeededNstrip", 15);
0284   preshowerSelection.add<double>("preshCalibPlaneX", 1.0);
0285   preshowerSelection.add<double>("preshCalibPlaneY", 0.7);
0286   preshowerSelection.add<double>("preshCalibGamma", 0.024);
0287   preshowerSelection.add<double>("preshCalibMIP", 9.0E-5);
0288   preshowerSelection.add<std::string>("debugLevelES", "");  // *** This is not needed and shoul be better removed !
0289   desc.add<edm::ParameterSetDescription>("preshowerSelection", preshowerSelection);
0290 
0291   desc.add<int>("debugLevel", 0);
0292   descriptions.add("hltRegionalEcalResonanceFilter", desc);
0293 }
0294 
0295 // ------------ method called to produce the data  ------------
0296 bool HLTRegionalEcalResonanceFilter::filter(edm::Event &iEvent, const edm::EventSetup &iSetup) {
0297   //Create empty output collections
0298   std::unique_ptr<EBRecHitCollection> selEBRecHitCollection(new EBRecHitCollection);
0299   //Create empty output collections
0300   std::unique_ptr<EERecHitCollection> selEERecHitCollection(new EERecHitCollection);
0301 
0302   ////all selected..
0303   vector<DetId> selectedEBDetIds;
0304   vector<DetId> selectedEEDetIds;
0305 
0306   auto const &pTopology = iSetup.getHandle(caloTopologyRecordToken_);
0307   const CaloSubdetectorTopology *topology_eb = pTopology->getSubdetectorTopology(DetId::Ecal, EcalBarrel);
0308   const CaloSubdetectorTopology *topology_ee = pTopology->getSubdetectorTopology(DetId::Ecal, EcalEndcap);
0309 
0310   auto const &geoHandle = iSetup.getHandle(caloGeometryRecordToken_);
0311   const CaloSubdetectorGeometry *geometry_es = geoHandle->getSubdetectorGeometry(DetId::Ecal, EcalPreshower);
0312 
0313   std::unique_ptr<CaloSubdetectorTopology> topology_es;
0314   if (geometry_es) {
0315     topology_es = std::make_unique<EcalPreshowerTopology>();
0316   } else {
0317     storeRecHitES_ = false;  ///if no preshower
0318   }
0319 
0320   ///get status from DB
0321   edm::ESHandle<EcalChannelStatus> csHandle;
0322   if (useDBStatus_)
0323     csHandle = iSetup.getHandle(ecalChannelStatusRcdToken_);
0324   const EcalChannelStatus &channelStatus = *csHandle;
0325 
0326   ///==============Start to process barrel part==================///
0327 
0328   Handle<EBRecHitCollection> barrelRecHitsHandle;
0329 
0330   iEvent.getByToken(barrelHitsToken_, barrelRecHitsHandle);
0331   if (!barrelRecHitsHandle.isValid()) {
0332     LogError("AlCaEcalResonanceProducer") << "AlCaEcalResonanceProducer: Error! can't get product barrel hits!";
0333   }
0334 
0335   const EcalRecHitCollection *hitCollection_eb = barrelRecHitsHandle.product();
0336 
0337   Handle<reco::BasicClusterCollection> barrelClustersHandle;
0338   iEvent.getByToken(barrelClustersToken_, barrelClustersHandle);
0339   if (!barrelClustersHandle.isValid()) {
0340     LogError("AlCaEcalResonanceProducer") << "AlCaEcalResonanceProducer: Error! can't get product barrel clusters!";
0341   }
0342   const reco::BasicClusterCollection *clusterCollection_eb = barrelClustersHandle.product();
0343   if (debug_ >= 1)
0344     LogDebug("AlCaEcalResonanceProducer")
0345         << " barrel_input_size:  run " << iEvent.id().run() << " event " << iEvent.id().event()
0346         << " nhitEB:  " << hitCollection_eb->size() << " nCluster: " << clusterCollection_eb->size();
0347 
0348   if (doSelBarrel_) {
0349     map<int, vector<EcalRecHit> > RecHits5x5_clus;  ///5x5 for selected pairs
0350     vector<int> indIsoClus;                         /// Iso cluster all , 5x5 rechit not yet done
0351     vector<int> indCandClus;                        ///good cluster all ,  5x5 rechit done already during the loop
0352     vector<int> indClusSelected;                    /// saved so far, all
0353 
0354     doSelection(EcalBarrel,
0355                 clusterCollection_eb,
0356                 hitCollection_eb,
0357                 channelStatus,
0358                 topology_eb,
0359                 RecHits5x5_clus,
0360                 indCandClus,
0361                 indIsoClus,
0362                 indClusSelected);
0363 
0364     ///Now save all rechits in the selected clusters
0365     for (int ind : indClusSelected) {
0366       auto it_bc3 = clusterCollection_eb->begin() + ind;
0367       const std::vector<std::pair<DetId, float> > &vid = it_bc3->hitsAndFractions();
0368       for (auto const &idItr : vid) {
0369         auto hit = hitCollection_eb->find(idItr.first);
0370         if (hit == hitCollection_eb->end())
0371           continue;  //this should not happen.
0372         selEBRecHitCollection->push_back(*hit);
0373         selectedEBDetIds.push_back(idItr.first);
0374       }
0375     }
0376 
0377     if (store5x5RecHitEB_) {
0378       ///stroe 5x5 of good clusters, 5x5 arleady got
0379       for (int ind : indCandClus) {
0380         vector<EcalRecHit> v5x5 = RecHits5x5_clus[ind];
0381         for (auto &n : v5x5) {
0382           DetId ed = n.id();
0383           auto itdet = find(selectedEBDetIds.begin(), selectedEBDetIds.end(), ed);
0384           if (itdet == selectedEBDetIds.end()) {
0385             selectedEBDetIds.push_back(ed);
0386             selEBRecHitCollection->push_back(n);
0387           }
0388         }
0389       }
0390 
0391       ///store 5x5 of Iso clusters, need to getWindow of 5x5
0392       for (int ind : indIsoClus) {
0393         auto it =
0394             find(indCandClus.begin(), indCandClus.end(), ind);  ///check if already saved in the good cluster vector
0395         if (it != indCandClus.end())
0396           continue;
0397 
0398         auto it_bc3 = clusterCollection_eb->begin() + ind;
0399         DetId seedId = it_bc3->seed();
0400         std::vector<DetId> clus_v5x5 = topology_eb->getWindow(seedId, 5, 5);
0401         for (std::vector<DetId>::const_iterator idItr = clus_v5x5.begin(); idItr != clus_v5x5.end(); idItr++) {
0402           DetId ed = *idItr;
0403           auto rit = hitCollection_eb->find(ed);
0404           if (rit == hitCollection_eb->end())
0405             continue;
0406           if (!checkStatusOfEcalRecHit(channelStatus, *rit))
0407             continue;
0408           auto itdet = find(selectedEBDetIds.begin(), selectedEBDetIds.end(), ed);
0409           if (itdet == selectedEBDetIds.end()) {
0410             selectedEBDetIds.push_back(ed);
0411             selEBRecHitCollection->push_back(*rit);
0412           }
0413         }
0414       }
0415     }
0416 
0417   }  // end of selection for eta/pi0->gg in the barrel
0418 
0419   int eb_collsize = selEBRecHitCollection->size();
0420 
0421   if (debug_ >= 1)
0422     LogDebug("AlCaEcalResonanceProducer") << " barrel_output_size_run " << iEvent.id().run() << " event "
0423                                           << iEvent.id().event() << " nEBSaved " << selEBRecHitCollection->size();
0424   ///==============End of  barrel ==================///
0425 
0426   //===============Start of Endcap =================/////
0427   ///get preshower rechits
0428   Handle<ESRecHitCollection> esRecHitsHandle;
0429   iEvent.getByToken(preshHitsToken_, esRecHitsHandle);
0430   if (!esRecHitsHandle.isValid()) {
0431     LogError("AlCaEcalResonanceProducer") << "AlCaEcalResonanceProducer: Error! can't get product esRecHit!";
0432   }
0433   const EcalRecHitCollection *hitCollection_es = esRecHitsHandle.product();
0434   // make a map of rechits:
0435   m_esrechit_map.clear();
0436   EcalRecHitCollection::const_iterator iter;
0437   for (iter = esRecHitsHandle->begin(); iter != esRecHitsHandle->end(); iter++) {
0438     //Make the map of DetID, EcalRecHit pairs
0439     m_esrechit_map.insert(std::make_pair(iter->id(), *iter));
0440   }
0441   // The set of used DetID's for a given event:
0442   m_used_strips.clear();
0443   std::unique_ptr<ESRecHitCollection> selESRecHitCollection(new ESRecHitCollection);
0444 
0445   Handle<EERecHitCollection> endcapRecHitsHandle;
0446   iEvent.getByToken(endcapHitsToken_, endcapRecHitsHandle);
0447   if (!endcapRecHitsHandle.isValid()) {
0448     LogError("AlCaEcalResonanceProducer") << "AlCaEcalResonanceProducer: Error! can't get product endcap hits!";
0449   }
0450   const EcalRecHitCollection *hitCollection_ee = endcapRecHitsHandle.product();
0451   Handle<reco::BasicClusterCollection> endcapClustersHandle;
0452   iEvent.getByToken(endcapClustersToken_, endcapClustersHandle);
0453   if (!endcapClustersHandle.isValid()) {
0454     LogError("AlCaEcalResonanceProducer") << "AlCaEcalResonanceProducer: Error! can't get product endcap clusters!";
0455   }
0456   const reco::BasicClusterCollection *clusterCollection_ee = endcapClustersHandle.product();
0457   if (debug_ >= 1)
0458     LogDebug("AlCaEcalResonanceProducer")
0459         << " endcap_input_size:  run " << iEvent.id().run() << " event " << iEvent.id().event()
0460         << " nhitEE:  " << hitCollection_ee->size() << " nhitES: " << hitCollection_es->size()
0461         << " nCluster: " << clusterCollection_ee->size();
0462 
0463   if (doSelEndcap_) {
0464     map<int, vector<EcalRecHit> > RecHits5x5_clus;  ///5x5 for selected pairs
0465     vector<int> indIsoClus;                         /// Iso cluster all , 5x5 rechit not yet done
0466     vector<int> indCandClus;                        ///good cluster all ,  5x5 rechit done already during the loop
0467     vector<int> indClusSelected;                    /// saved so far, all
0468 
0469     doSelection(EcalEndcap,
0470                 clusterCollection_ee,
0471                 hitCollection_ee,
0472                 channelStatus,
0473                 topology_ee,
0474                 RecHits5x5_clus,
0475                 indCandClus,
0476                 indIsoClus,
0477                 indClusSelected);
0478 
0479     ///Now save all rechits in the selected clusters
0480     for (int ind : indClusSelected) {
0481       auto it_bc3 = clusterCollection_ee->begin() + ind;
0482       const std::vector<std::pair<DetId, float> > &vid = it_bc3->hitsAndFractions();
0483       for (auto const &idItr : vid) {
0484         auto hit = hitCollection_ee->find(idItr.first);
0485         if (hit == hitCollection_ee->end())
0486           continue;  //this should not happen.
0487         selEERecHitCollection->push_back(*hit);
0488         selectedEEDetIds.push_back(idItr.first);
0489       }
0490       ///save preshower rechits
0491       if (storeRecHitES_) {
0492         std::set<DetId> used_strips_before = m_used_strips;
0493         makeClusterES(it_bc3->x(), it_bc3->y(), it_bc3->z(), geometry_es, topology_es.get());
0494         auto ites = m_used_strips.begin();
0495         for (; ites != m_used_strips.end(); ++ites) {
0496           ESDetId d1 = ESDetId(*ites);
0497           auto ites2 = find(used_strips_before.begin(), used_strips_before.end(), d1);
0498           if ((ites2 == used_strips_before.end())) {
0499             auto itmap = m_esrechit_map.find(d1);
0500             selESRecHitCollection->push_back(itmap->second);
0501           }
0502         }
0503       }
0504     }
0505 
0506     if (store5x5RecHitEE_) {
0507       ///store 5x5 of good clusters, 5x5 arleady got
0508       for (int ind : indCandClus) {
0509         vector<EcalRecHit> v5x5 = RecHits5x5_clus[ind];
0510         for (auto &n : v5x5) {
0511           DetId ed = n.id();
0512           auto itdet = find(selectedEEDetIds.begin(), selectedEEDetIds.end(), ed);
0513           if (itdet == selectedEEDetIds.end()) {
0514             selectedEEDetIds.push_back(ed);
0515             selEERecHitCollection->push_back(n);
0516           }
0517         }
0518       }
0519 
0520       ///store 5x5 of Iso clusters, need to getWindow of 5x5
0521       for (int ind : indIsoClus) {
0522         auto it =
0523             find(indCandClus.begin(), indCandClus.end(), ind);  ///check if already saved in the good cluster vector
0524         if (it != indCandClus.end())
0525           continue;
0526 
0527         auto it_bc3 = clusterCollection_ee->begin() + ind;
0528         DetId seedId = it_bc3->seed();
0529         std::vector<DetId> clus_v5x5 = topology_ee->getWindow(seedId, 5, 5);
0530         for (std::vector<DetId>::const_iterator idItr = clus_v5x5.begin(); idItr != clus_v5x5.end(); idItr++) {
0531           DetId ed = *idItr;
0532           auto rit = hitCollection_ee->find(ed);
0533           if (rit == hitCollection_ee->end())
0534             continue;
0535           if (!checkStatusOfEcalRecHit(channelStatus, *rit))
0536             continue;
0537           auto itdet = find(selectedEEDetIds.begin(), selectedEEDetIds.end(), ed);
0538           if (itdet == selectedEEDetIds.end()) {
0539             selectedEEDetIds.push_back(ed);
0540             selEERecHitCollection->push_back(*rit);
0541           }
0542         }
0543       }
0544     }
0545 
0546   }  // end of selection for eta/pi0->gg in the endcap
0547 
0548   ////==============End of endcap ===============///
0549 
0550   if (debug_ >= 1)
0551     LogDebug("AlCaEcalResonanceProducer")
0552         << " endcap_output_size run " << iEvent.id().run() << "_" << iEvent.id().event() << " nEESaved "
0553         << selEERecHitCollection->size() << " nESSaved: " << selESRecHitCollection->size();
0554 
0555   //Put selected information in the event
0556   int ee_collsize = selEERecHitCollection->size();
0557 
0558   ////Now put into events selected rechits.
0559   if (doSelBarrel_) {
0560     iEvent.put(std::move(selEBRecHitCollection), BarrelHits_);
0561   }
0562   if (doSelEndcap_) {
0563     iEvent.put(std::move(selEERecHitCollection), EndcapHits_);
0564     if (storeRecHitES_) {
0565       iEvent.put(std::move(selESRecHitCollection), ESHits_);
0566     }
0567   }
0568 
0569   return (eb_collsize > 1 or ee_collsize > 1);
0570 }
0571 
0572 void HLTRegionalEcalResonanceFilter::doSelection(
0573     int detector,
0574     const reco::BasicClusterCollection *clusterCollection,
0575     const EcalRecHitCollection *hitCollection,
0576     const EcalChannelStatus &channelStatus,
0577     const CaloSubdetectorTopology *topology_p,
0578     map<int, vector<EcalRecHit> > &RecHits5x5_clus,
0579     vector<int> &indCandClus,     ///good cluster all ,  5x5 rechit done already during the loop
0580     vector<int> &indIsoClus,      /// Iso cluster all , 5x5 rechit not yet done
0581     vector<int> &indClusSelected  /// saved so far, all
0582 ) {
0583   vector<int> indClusPi0Candidates;  ///those clusters identified as pi0s
0584   if (detector == EcalBarrel && removePi0CandidatesForEta_) {
0585     for (auto it_bc = clusterCollection->begin(); it_bc != clusterCollection->end(); it_bc++) {
0586       for (auto it_bc2 = it_bc + 1; it_bc2 != clusterCollection->end(); it_bc2++) {
0587         float m_pair, pt_pair, eta_pair, phi_pair;
0588         calcPaircluster(*it_bc, *it_bc2, m_pair, pt_pair, eta_pair, phi_pair);
0589         if (m_pair > massLowPi0Cand_ && m_pair < massHighPi0Cand_) {
0590           int indtmp[2] = {int(it_bc - clusterCollection->begin()), int(it_bc2 - clusterCollection->begin())};
0591           for (int &k : indtmp) {
0592             auto it = find(indClusPi0Candidates.begin(), indClusPi0Candidates.end(), k);
0593             if (it == indClusPi0Candidates.end())
0594               indClusPi0Candidates.push_back(k);
0595           }
0596         }
0597       }
0598     }  //end of loop over finding pi0's clusters
0599   }
0600 
0601   //mass cut
0602   double m_minCut = seleMinvMinBarrel_;
0603   double m_maxCut = seleMinvMaxBarrel_;
0604 
0605   //isolation
0606   double ptminforIso = ptMinForIsolation_;
0607   double isoBeltdrCut = seleBeltDR_;
0608   double isoBeltdetaCut = seleBeltDeta_;
0609 
0610   //other
0611   double s9s25Cut = seleS9S25Gamma_;
0612   bool store5x5 = store5x5RecHitEB_;
0613 
0614   //double isoCut = seleIso_;  //old non-region filter
0615   //double s4s9Cut = seleS4S9Gamma_;  // old non-region filter
0616 
0617   if (detector == EcalEndcap) {
0618     //mass cuts
0619     m_minCut = seleMinvMinEndCap_;
0620     m_maxCut = seleMinvMaxEndCap_;
0621 
0622     //isolation
0623     ptminforIso = ptMinForIsolationEndCap_;
0624     isoBeltdrCut = seleBeltDREndCap_;
0625     isoBeltdetaCut = seleBeltDetaEndCap_;
0626 
0627     //other
0628     s9s25Cut = seleS9S25GammaEndCap_;
0629     store5x5 = store5x5RecHitEE_;
0630 
0631     //isoCut = seleIsoEndCap_;  // old non-region filter
0632     //s4s9Cut = seleS4S9GammaEndCap_;  //old non-region filter
0633   }
0634 
0635   map<int, bool>
0636       passShowerShape_clus;  //if this cluster passed showershape cut, so no need to compute the quantity again for each loop
0637 
0638   for (auto it_bc = clusterCollection->begin(); it_bc != clusterCollection->end(); it_bc++) {
0639     if (detector == EcalBarrel && removePi0CandidatesForEta_) {
0640       auto it = find(indClusPi0Candidates.begin(), indClusPi0Candidates.end(), int(it_bc - clusterCollection->begin()));
0641       if (it != indClusPi0Candidates.end())
0642         continue;
0643     }
0644 
0645     float en1 = it_bc->energy();
0646     math::XYZVector v1(it_bc->position());  // set vector as cluster position
0647     v1 *=
0648         (en1 /
0649          v1.R());  // rescale vector's magnitude to the energy in order to get momentum vector (assuming massless particles)
0650     float pt1 = v1.Rho();  // Rho is equivalent to Pt when using XYZVector
0651 
0652     int ind1 = int(it_bc - clusterCollection->begin());
0653     auto itmap = passShowerShape_clus.find(ind1);
0654     if (itmap != passShowerShape_clus.end()) {
0655       if (itmap->second == false) {
0656         continue;
0657       }
0658     }
0659 
0660     for (auto it_bc2 = it_bc + 1; it_bc2 != clusterCollection->end(); it_bc2++) {
0661       if (detector == EcalBarrel && removePi0CandidatesForEta_) {
0662         auto it =
0663             find(indClusPi0Candidates.begin(), indClusPi0Candidates.end(), int(it_bc2 - clusterCollection->begin()));
0664         if (it != indClusPi0Candidates.end())
0665           continue;
0666       }
0667 
0668       float en2 = it_bc2->energy();
0669       math::XYZVector v2(it_bc2->position());  // set vector as cluster position
0670       v2 *=
0671           (en2 /
0672            v2.R());  // rescale vector's magnitude to the energy in order to get momentum vector (assuming massless particles)
0673       float pt2 = v2.Rho();  // Rho is equivalent to Pt when using XYZVector
0674 
0675       int ind2 = int(it_bc2 - clusterCollection->begin());
0676       auto itmap = passShowerShape_clus.find(ind2);
0677       if (itmap != passShowerShape_clus.end()) {
0678         if (itmap->second == false) {
0679           continue;
0680         }
0681       }
0682 
0683       float m_pair, pt_pair, eta_pair, phi_pair;
0684       calcPaircluster(*it_bc, *it_bc2, m_pair, pt_pair, eta_pair, phi_pair);
0685 
0686       /// pt & Pt pair Cut
0687       float ptmin = pt1 < pt2 ? pt1 : pt2;
0688       float etapair = fabs(eta_pair);
0689 
0690       //-------------------------------------
0691       // Region Based Kinematic Cuts: pt of the diphoton system and pt of each photon
0692       //-------------------------------------
0693       if (detector == EcalBarrel) {  // BARREL
0694         //EB region 1
0695         if (etapair <= region1_Barrel_) {  //EB region 1
0696           if (ptmin < selePtGammaBarrel_region1_ || pt_pair < selePtPairBarrel_region1_)
0697             continue;
0698         }
0699         //EB region 2
0700         else {
0701           if (ptmin < selePtGammaBarrel_region2_ || pt_pair < selePtPairBarrel_region2_)
0702             continue;
0703         }
0704       } else {  // ENDCAP
0705         //EE region 1
0706         if (etapair <= region1_EndCap_) {
0707           if (ptmin < selePtGammaEndCap_region1_ || pt_pair < selePtPairEndCap_region1_)
0708             continue;
0709         }
0710         //EE region 2
0711         else if (etapair <= region2_EndCap_) {
0712           if (ptmin < selePtGammaEndCap_region2_ || pt_pair < selePtPairEndCap_region2_)
0713             continue;
0714         }
0715         //EE region 3
0716         else {
0717           if (ptmin < selePtGammaEndCap_region3_ || pt_pair < selePtPairEndCap_region3_)
0718             continue;
0719           if (pt_pair > selePtPairMaxEndCap_region3_)
0720             continue;  // there is also a possible maximum pt for the very forward region
0721         }
0722       }
0723 
0724       /// Mass window Cut
0725       if (m_pair < m_minCut || m_pair > m_maxCut)
0726         continue;
0727 
0728       //// Loop on cluster to measure isolation:
0729       vector<int> IsoClus;
0730       IsoClus.push_back(ind1);  //first two are good clusters
0731       IsoClus.push_back(ind2);
0732 
0733       float Iso = 0;
0734       for (auto it_bc3 = clusterCollection->begin(); it_bc3 != clusterCollection->end(); it_bc3++) {
0735         if (it_bc3->seed() == it_bc->seed() || it_bc3->seed() == it_bc2->seed())
0736           continue;
0737         float drcl = std::hypot(eta_pair - it_bc3->eta(), deltaPhi(phi_pair, it_bc3->phi()));
0738         float dretacl = fabs(eta_pair - it_bc3->eta());
0739         if (drcl > isoBeltdrCut || dretacl > isoBeltdetaCut)
0740           continue;
0741         float pt3 = it_bc3->energy() * sin(it_bc3->position().theta());
0742         if (pt3 < ptminforIso)
0743           continue;
0744         Iso += pt3;
0745         int ind3 = int(it_bc3 - clusterCollection->begin());  /// remember which Iso cluster used
0746         IsoClus.push_back(ind3);
0747       }
0748 
0749       //-------------------------------------
0750       // Region Based Nxtal Cut: Nxtal of the clusters
0751       //-------------------------------------
0752 
0753       uint32_t Nxtal_pho1 = it_bc->size();
0754       uint32_t Nxtal_pho2 = it_bc2->size();
0755 
0756       // BARREL
0757       if (detector == EcalBarrel) {
0758         //EB region 1
0759         if (etapair <= region1_Barrel_) {  //EB region 1
0760           if (Nxtal_pho1 < seleNxtalBarrel_region1_ || Nxtal_pho2 < seleNxtalBarrel_region1_)
0761             continue;
0762         }
0763         //EB region 2
0764         else {
0765           if (Nxtal_pho1 < seleNxtalBarrel_region2_ || Nxtal_pho2 < seleNxtalBarrel_region2_)
0766             continue;
0767         }
0768       }
0769       // ENDCAP
0770       else {
0771         //EE region 1
0772         if (etapair <= region1_EndCap_) {
0773           if (Nxtal_pho1 < seleNxtalEndCap_region1_ || Nxtal_pho2 < seleNxtalEndCap_region1_)
0774             continue;
0775         }
0776         //EE region 2
0777         else if (etapair <= region2_EndCap_) {
0778           if (Nxtal_pho1 < seleNxtalEndCap_region2_ || Nxtal_pho2 < seleNxtalEndCap_region2_)
0779             continue;
0780         }
0781         //EE region 3
0782         else {
0783           if (Nxtal_pho1 < seleNxtalEndCap_region3_ || Nxtal_pho2 < seleNxtalEndCap_region3_)
0784             continue;
0785         }
0786       }
0787 
0788       //-------------------------------------
0789       // Region Based Isolation Cut: pt of the diphoton system and pt of each photon
0790       //-------------------------------------
0791       float iso_val = Iso / pt_pair;
0792       // BARREL
0793       if (detector == EcalBarrel) {
0794         //EB region 1
0795         if (etapair <= region1_Barrel_) {  //EB region 1
0796           if (iso_val > seleIsoBarrel_region1_)
0797             continue;
0798         }
0799         //EB region 2
0800         else {
0801           if (iso_val > seleIsoBarrel_region2_)
0802             continue;
0803         }
0804       }
0805       // ENDCAP
0806       else {
0807         //EE region 1
0808         if (etapair <= region1_EndCap_) {
0809           if (iso_val > seleIsoEndCap_region1_)
0810             continue;
0811         }
0812         //EE region 2
0813         else if (etapair <= region2_EndCap_) {
0814           if (iso_val > seleIsoEndCap_region2_)
0815             continue;
0816         }
0817         //EE region 3
0818         else {
0819           if (iso_val > seleIsoEndCap_region3_)
0820             continue;
0821         }
0822       }
0823       //-------------------------------------
0824       //if(Iso/pt_pair > isoCut) continue;  //old non-regional cut
0825 
0826       //-------------------------------------
0827       // Region based ShowerShape Cut: S4S9 Cut with possible S9S25 configuration
0828       //-------------------------------------
0829 
0830       bool failShowerShape = false;
0831       for (int n = 0; n < 2; n++) {
0832         int ind = IsoClus[n];
0833         auto it_bc3 = clusterCollection->begin() + ind;
0834         auto itmap = passShowerShape_clus.find(ind);
0835 
0836         if (itmap != passShowerShape_clus.end()) {  // if we havent reached the end
0837           if (itmap->second == false) {
0838             failShowerShape = true;
0839             n = 2;  //exit the loop
0840           }
0841         } else {
0842           vector<EcalRecHit> RecHitsCluster_5x5;
0843           float res[3];
0844 
0845           //build the shower shape variables
0846           calcShowerShape(*it_bc3, channelStatus, hitCollection, topology_p, store5x5, RecHitsCluster_5x5, res);
0847           float s4s9 = res[1] > 0 ? res[0] / res[1] : -1;
0848           float s9s25 = res[2] > 0 ? res[1] / res[2] : -1;
0849 
0850           bool passed = false;
0851 
0852           // BARREL
0853           if (detector == EcalBarrel) {
0854             //EB region 1
0855             if (etapair <= region1_Barrel_) {
0856               passed = s4s9 > seleS4S9GammaBarrel_region1_;
0857             }
0858             //EB region 2
0859             else {
0860               passed = s4s9 > seleS4S9GammaBarrel_region2_;
0861             }
0862           }
0863           // ENDCAP
0864           else {
0865             //EE region 1
0866             if (etapair <= region1_EndCap_) {
0867               passed = s4s9 > seleS4S9GammaEndCap_region1_;
0868             }
0869             //EE region 2
0870             else if (etapair <= region2_EndCap_) {
0871               passed = s4s9 > seleS4S9GammaEndCap_region2_;
0872             }
0873             //EE region 3
0874             else {
0875               passed = s4s9 > seleS4S9GammaEndCap_region3_;
0876             }
0877           }
0878 
0879           //apply the S9S25 cut as well
0880           passed = passed && (s9s25 > s9s25Cut);
0881           passShowerShape_clus.insert(pair<int, bool>(ind, passed));
0882 
0883           if (!passed) {
0884             failShowerShape = true;
0885             n = 2;  //exit the loop
0886           } else {
0887             RecHits5x5_clus.insert(pair<int, vector<EcalRecHit> >(ind, RecHitsCluster_5x5));
0888           }
0889         }
0890       }
0891 
0892       if (failShowerShape == true)
0893         continue;  //if any of the two clusters fail shower shape
0894 
0895       ///Save two good clusters' index( plus Iso cluster )  if not yet saved
0896       for (unsigned int iii = 0; iii < IsoClus.size(); iii++) {
0897         int ind = IsoClus[iii];
0898 
0899         if (iii < 2) {
0900           auto it = find(indCandClus.begin(), indCandClus.end(), ind);  ///good cluster
0901           if (it == indCandClus.end())
0902             indCandClus.push_back(ind);
0903           else
0904             continue;
0905         } else {
0906           auto it = find(indIsoClus.begin(), indIsoClus.end(), ind);  //iso cluster
0907           if (it == indIsoClus.end())
0908             indIsoClus.push_back(ind);
0909           else
0910             continue;
0911         }
0912 
0913         auto it = find(indClusSelected.begin(), indClusSelected.end(), ind);
0914         if (it != indClusSelected.end())
0915           continue;
0916         indClusSelected.push_back(ind);
0917       }
0918     }
0919   }
0920 }
0921 
0922 void HLTRegionalEcalResonanceFilter::convxtalid(Int_t &nphi, Int_t &neta) {
0923   // Barrel only
0924   // Output nphi 0...359; neta 0...84; nside=+1 (for eta>0), or 0 (for eta<0).
0925   // neta will be [-85,-1] , or [0,84], the minus sign indicates the z<0 side.
0926 
0927   if (neta > 0)
0928     neta -= 1;
0929   if (nphi > 359)
0930     nphi = nphi - 360;
0931 
0932   // final check
0933   if (nphi > 359 || nphi < 0 || neta < -85 || neta > 84) {
0934     LogError("") << " unexpected fatal error in HLTEcalResonanceFilter::convxtalid " << nphi << " " << neta << " "
0935                  << std::endl;
0936   }
0937 }  //end of convxtalid
0938 
0939 int HLTRegionalEcalResonanceFilter::diff_neta_s(Int_t neta1, Int_t neta2) {
0940   Int_t mdiff;
0941   mdiff = (neta1 - neta2);
0942   return mdiff;
0943 }
0944 
0945 // Calculate the distance in xtals taking into account the periodicity of the Barrel
0946 int HLTRegionalEcalResonanceFilter::diff_nphi_s(Int_t nphi1, Int_t nphi2) {
0947   Int_t mdiff;
0948   if (std::abs(nphi1 - nphi2) < (360 - std::abs(nphi1 - nphi2))) {
0949     mdiff = nphi1 - nphi2;
0950   } else {
0951     mdiff = 360 - std::abs(nphi1 - nphi2);
0952     if (nphi1 > nphi2)
0953       mdiff = -mdiff;
0954   }
0955   return mdiff;
0956 }
0957 
0958 void HLTRegionalEcalResonanceFilter::calcShowerShape(const reco::BasicCluster &bc,
0959                                                      const EcalChannelStatus &channelStatus,
0960                                                      const EcalRecHitCollection *recHits,
0961                                                      const CaloSubdetectorTopology *topology_p,
0962                                                      bool calc5x5,
0963                                                      vector<EcalRecHit> &rechit5x5,
0964                                                      float res[]) {
0965   const std::vector<std::pair<DetId, float> > &vid = bc.hitsAndFractions();
0966 
0967   float e2x2[4] = {0};
0968   float e3x3 = 0;
0969   float e5x5 = 0;
0970   int seedx;
0971   int seedy;
0972   DetId seedId = bc.seed();
0973 
0974   bool InBarrel = true;
0975   if (seedId.subdetId() == EcalBarrel) {
0976     EBDetId ebd = EBDetId(seedId);
0977     seedx = ebd.ieta();
0978     seedy = ebd.iphi();
0979     convxtalid(seedy, seedx);
0980   } else {
0981     InBarrel = false;
0982     EEDetId eed = EEDetId(seedId);
0983     seedx = eed.ix();
0984     seedy = eed.iy();
0985   }
0986   int x, y, dx, dy;
0987 
0988   if (calc5x5) {
0989     rechit5x5.clear();
0990     std::vector<DetId> clus_v5x5;
0991     clus_v5x5 = topology_p->getWindow(seedId, 5, 5);
0992 
0993     for (std::vector<DetId>::const_iterator idItr = clus_v5x5.begin(); idItr != clus_v5x5.end(); idItr++) {
0994       DetId ed = *idItr;
0995       if (InBarrel == true) {
0996         EBDetId ebd = EBDetId(ed);
0997         x = ebd.ieta();
0998         y = ebd.iphi();
0999         convxtalid(y, x);
1000         dx = diff_neta_s(x, seedx);
1001         dy = diff_nphi_s(y, seedy);
1002       } else {
1003         EEDetId eed = EEDetId(ed);
1004         x = eed.ix();
1005         y = eed.iy();
1006         dx = x - seedx;
1007         dy = y - seedy;
1008       }
1009       auto rit = recHits->find(ed);
1010       if (rit == recHits->end())
1011         continue;
1012       if (!checkStatusOfEcalRecHit(channelStatus, *rit))
1013         continue;
1014 
1015       float energy = (*rit).energy();
1016       e5x5 += energy;
1017 
1018       auto idItrF = std::find(vid.begin(), vid.end(), std::make_pair(ed, 1.0f));  ///has to add "f", make it float
1019       if (idItrF == vid.end()) {  ///only store those not belonging to this cluster
1020         rechit5x5.push_back(*rit);
1021       } else {  /// S4, S9 are defined inside the cluster, the same as below.
1022         if (std::abs(dx) <= 1 && std::abs(dy) <= 1) {
1023           if (dx <= 0 && dy <= 0)
1024             e2x2[0] += energy;
1025           if (dx >= 0 && dy <= 0)
1026             e2x2[1] += energy;
1027           if (dx <= 0 && dy >= 0)
1028             e2x2[2] += energy;
1029           if (dx >= 0 && dy >= 0)
1030             e2x2[3] += energy;
1031           e3x3 += energy;
1032         }
1033       }
1034     }
1035 
1036   } else {
1037     for (auto const &idItr : vid) {
1038       DetId ed = idItr.first;
1039       if (InBarrel == true) {
1040         EBDetId ebd = EBDetId(ed);
1041         x = ebd.ieta();
1042         y = ebd.iphi();
1043         convxtalid(y, x);
1044         dx = diff_neta_s(x, seedx);
1045         dy = diff_nphi_s(y, seedy);
1046       } else {
1047         EEDetId eed = EEDetId(ed);
1048         x = eed.ix();
1049         y = eed.iy();
1050         dx = x - seedx;
1051         dy = y - seedy;
1052       }
1053       auto rit = recHits->find(ed);
1054       if (rit == recHits->end()) {
1055         continue;
1056       }
1057 
1058       float energy = (*rit).energy();
1059       if (std::abs(dx) <= 1 && std::abs(dy) <= 1) {
1060         if (dx <= 0 && dy <= 0)
1061           e2x2[0] += energy;
1062         if (dx >= 0 && dy <= 0)
1063           e2x2[1] += energy;
1064         if (dx <= 0 && dy >= 0)
1065           e2x2[2] += energy;
1066         if (dx >= 0 && dy >= 0)
1067           e2x2[3] += energy;
1068         e3x3 += energy;
1069       }
1070     }
1071     e5x5 = e3x3;  ///if not asked to calculte 5x5, then just make e5x5 = e3x3
1072   }
1073 
1074   ///e2x2
1075   res[0] = *max_element(e2x2, e2x2 + 4);
1076   res[1] = e3x3;
1077   res[2] = e5x5;
1078 }
1079 
1080 void HLTRegionalEcalResonanceFilter::calcPaircluster(const reco::BasicCluster &bc1,
1081                                                      const reco::BasicCluster &bc2,
1082                                                      float &m_pair,
1083                                                      float &pt_pair,
1084                                                      float &eta_pair,
1085                                                      float &phi_pair) {
1086   // use XYZVector instead of TLorentzVector to make things faster (and initialize with cartesian coordinates).
1087   // We are interested in the momentum vector:  however, we start from cartesian coordinates to get the vector direction,
1088   // then we set the vector's magnitude to obtain momentum coordinates. The magnitude we set is equal to the particle's energy.
1089   // We can do this because, assuming massless particles (or negligible mass), the magnitude of the momentum vector is given by the energy.
1090 
1091   math::XYZVector v1(bc1.position());
1092   float en1 = bc1.energy();
1093   float scaleFactor = en1 / v1.R();  // XYZVector::R() returns sqrt(Mag2()), where Mag2()= fx*fx + fy*fy + fz*fz
1094   // here I'm assuming that the vector initial magnitude is always different from 0 (the cluster must be somewhere in the detector, so the distance is greater than 0)
1095   v1 *= scaleFactor;
1096 
1097   // vsum would be v1 + v2, but instead of declaring both v2 and vsum, just declare vsum, initialize as if it is v2 and then sum v1.
1098   math::XYZVector vsum(bc2.position());
1099   // define energy sum initializing it to energy2, so that we can use it before summing energy1
1100   float energysum = bc2.energy();
1101   scaleFactor = energysum / vsum.R();
1102   vsum *= scaleFactor;
1103 
1104   vsum += v1;
1105   // now sum the energy of the second basic cluster to get total energy
1106   energysum += en1;
1107 
1108   // finally, assign values
1109   m_pair = sqrt(energysum * energysum - vsum.Mag2());  // M_pi0 = sqrt(E_pi0^2 - |p_pi0|^2)
1110   pt_pair = vsum.Rho();  // Rho method is the equivalent of Pt: returns sqrt( fx*fx + fy*fy )
1111   eta_pair = vsum.Eta();
1112   phi_pair = vsum.Phi();
1113 }
1114 
1115 bool HLTRegionalEcalResonanceFilter::checkStatusOfEcalRecHit(const EcalChannelStatus &channelStatus,
1116                                                              const EcalRecHit &rh) {
1117   if (useRecoFlag_) {  ///from recoFlag()
1118     int flag = rh.recoFlag();
1119     if (flagLevelRecHitsToUse_ == 0) {  ///good
1120       if (flag != 0)
1121         return false;
1122     } else if (flagLevelRecHitsToUse_ == 1) {  ///good || PoorCalib
1123       if (flag != 0 && flag != 4)
1124         return false;
1125     } else if (flagLevelRecHitsToUse_ == 2) {  ///good || PoorCalib || LeadingEdgeRecovered || kNeighboursRecovered,
1126       if (flag != 0 && flag != 4 && flag != 6 && flag != 7)
1127         return false;
1128     }
1129   }
1130   if (useDBStatus_) {  //// from DB
1131     int status = int(channelStatus[rh.id().rawId()].getStatusCode());
1132     if (status > statusLevelRecHitsToUse_)
1133       return false;
1134   }
1135 
1136   return true;
1137 }
1138 
1139 void HLTRegionalEcalResonanceFilter::makeClusterES(
1140     float x, float y, float z, const CaloSubdetectorGeometry *geometry_es, const CaloSubdetectorTopology *topology_es) {
1141   ///get assosicated ES clusters of this endcap cluster
1142   const GlobalPoint point(x, y, z);
1143   DetId tmp1 = (dynamic_cast<const EcalPreshowerGeometry *>(geometry_es))->getClosestCellInPlane(point, 1);
1144   DetId tmp2 = (dynamic_cast<const EcalPreshowerGeometry *>(geometry_es))->getClosestCellInPlane(point, 2);
1145   ESDetId strip1 = (tmp1 == DetId(0)) ? ESDetId(0) : ESDetId(tmp1);
1146   ESDetId strip2 = (tmp2 == DetId(0)) ? ESDetId(0) : ESDetId(tmp2);
1147 
1148   // Get ES clusters (found by the PreshSeeded algorithm) associated with a given EE cluster.
1149   for (int i2 = 0; i2 < preshNclust_; i2++) {
1150     reco::PreshowerCluster cl1 =
1151         presh_algo_->makeOneCluster(strip1, &m_used_strips, &m_esrechit_map, geometry_es, topology_es);
1152     reco::PreshowerCluster cl2 =
1153         presh_algo_->makeOneCluster(strip2, &m_used_strips, &m_esrechit_map, geometry_es, topology_es);
1154   }  // end of cycle over ES clusters
1155 }
1156 
1157 // declare this class as a framework plugin
1158 #include "FWCore/Framework/interface/MakerMacros.h"
1159 DEFINE_FWK_MODULE(HLTRegionalEcalResonanceFilter);
[ 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!