Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:18: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   if (eb_collsize < 2 && ee_collsize < 2)
0559     return false;
0560 
0561   ////Now put into events selected rechits.
0562   if (doSelBarrel_) {
0563     iEvent.put(std::move(selEBRecHitCollection), BarrelHits_);
0564   }
0565   if (doSelEndcap_) {
0566     iEvent.put(std::move(selEERecHitCollection), EndcapHits_);
0567     if (storeRecHitES_) {
0568       iEvent.put(std::move(selESRecHitCollection), ESHits_);
0569     }
0570   }
0571 
0572   return true;
0573 }
0574 
0575 void HLTRegionalEcalResonanceFilter::doSelection(
0576     int detector,
0577     const reco::BasicClusterCollection *clusterCollection,
0578     const EcalRecHitCollection *hitCollection,
0579     const EcalChannelStatus &channelStatus,
0580     const CaloSubdetectorTopology *topology_p,
0581     map<int, vector<EcalRecHit> > &RecHits5x5_clus,
0582     vector<int> &indCandClus,     ///good cluster all ,  5x5 rechit done already during the loop
0583     vector<int> &indIsoClus,      /// Iso cluster all , 5x5 rechit not yet done
0584     vector<int> &indClusSelected  /// saved so far, all
0585 ) {
0586   vector<int> indClusPi0Candidates;  ///those clusters identified as pi0s
0587   if (detector == EcalBarrel && removePi0CandidatesForEta_) {
0588     for (auto it_bc = clusterCollection->begin(); it_bc != clusterCollection->end(); it_bc++) {
0589       for (auto it_bc2 = it_bc + 1; it_bc2 != clusterCollection->end(); it_bc2++) {
0590         float m_pair, pt_pair, eta_pair, phi_pair;
0591         calcPaircluster(*it_bc, *it_bc2, m_pair, pt_pair, eta_pair, phi_pair);
0592         if (m_pair > massLowPi0Cand_ && m_pair < massHighPi0Cand_) {
0593           int indtmp[2] = {int(it_bc - clusterCollection->begin()), int(it_bc2 - clusterCollection->begin())};
0594           for (int &k : indtmp) {
0595             auto it = find(indClusPi0Candidates.begin(), indClusPi0Candidates.end(), k);
0596             if (it == indClusPi0Candidates.end())
0597               indClusPi0Candidates.push_back(k);
0598           }
0599         }
0600       }
0601     }  //end of loop over finding pi0's clusters
0602   }
0603 
0604   //mass cut
0605   double m_minCut = seleMinvMinBarrel_;
0606   double m_maxCut = seleMinvMaxBarrel_;
0607 
0608   //isolation
0609   double ptminforIso = ptMinForIsolation_;
0610   double isoBeltdrCut = seleBeltDR_;
0611   double isoBeltdetaCut = seleBeltDeta_;
0612 
0613   //other
0614   double s9s25Cut = seleS9S25Gamma_;
0615   bool store5x5 = store5x5RecHitEB_;
0616 
0617   //double isoCut = seleIso_;  //old non-region filter
0618   //double s4s9Cut = seleS4S9Gamma_;  // old non-region filter
0619 
0620   if (detector == EcalEndcap) {
0621     //mass cuts
0622     m_minCut = seleMinvMinEndCap_;
0623     m_maxCut = seleMinvMaxEndCap_;
0624 
0625     //isolation
0626     ptminforIso = ptMinForIsolationEndCap_;
0627     isoBeltdrCut = seleBeltDREndCap_;
0628     isoBeltdetaCut = seleBeltDetaEndCap_;
0629 
0630     //other
0631     s9s25Cut = seleS9S25GammaEndCap_;
0632     store5x5 = store5x5RecHitEE_;
0633 
0634     //isoCut = seleIsoEndCap_;  // old non-region filter
0635     //s4s9Cut = seleS4S9GammaEndCap_;  //old non-region filter
0636   }
0637 
0638   map<int, bool>
0639       passShowerShape_clus;  //if this cluster passed showershape cut, so no need to compute the quantity again for each loop
0640 
0641   for (auto it_bc = clusterCollection->begin(); it_bc != clusterCollection->end(); it_bc++) {
0642     if (detector == EcalBarrel && removePi0CandidatesForEta_) {
0643       auto it = find(indClusPi0Candidates.begin(), indClusPi0Candidates.end(), int(it_bc - clusterCollection->begin()));
0644       if (it != indClusPi0Candidates.end())
0645         continue;
0646     }
0647 
0648     float en1 = it_bc->energy();
0649     math::XYZVector v1(it_bc->position());  // set vector as cluster position
0650     v1 *=
0651         (en1 /
0652          v1.R());  // rescale vector's magnitude to the energy in order to get momentum vector (assuming massless particles)
0653     float pt1 = v1.Rho();  // Rho is equivalent to Pt when using XYZVector
0654 
0655     int ind1 = int(it_bc - clusterCollection->begin());
0656     auto itmap = passShowerShape_clus.find(ind1);
0657     if (itmap != passShowerShape_clus.end()) {
0658       if (itmap->second == false) {
0659         continue;
0660       }
0661     }
0662 
0663     for (auto it_bc2 = it_bc + 1; it_bc2 != clusterCollection->end(); it_bc2++) {
0664       if (detector == EcalBarrel && removePi0CandidatesForEta_) {
0665         auto it =
0666             find(indClusPi0Candidates.begin(), indClusPi0Candidates.end(), int(it_bc2 - clusterCollection->begin()));
0667         if (it != indClusPi0Candidates.end())
0668           continue;
0669       }
0670 
0671       float en2 = it_bc2->energy();
0672       math::XYZVector v2(it_bc2->position());  // set vector as cluster position
0673       v2 *=
0674           (en2 /
0675            v2.R());  // rescale vector's magnitude to the energy in order to get momentum vector (assuming massless particles)
0676       float pt2 = v2.Rho();  // Rho is equivalent to Pt when using XYZVector
0677 
0678       int ind2 = int(it_bc2 - clusterCollection->begin());
0679       auto itmap = passShowerShape_clus.find(ind2);
0680       if (itmap != passShowerShape_clus.end()) {
0681         if (itmap->second == false) {
0682           continue;
0683         }
0684       }
0685 
0686       float m_pair, pt_pair, eta_pair, phi_pair;
0687       calcPaircluster(*it_bc, *it_bc2, m_pair, pt_pair, eta_pair, phi_pair);
0688 
0689       /// pt & Pt pair Cut
0690       float ptmin = pt1 < pt2 ? pt1 : pt2;
0691       float etapair = fabs(eta_pair);
0692 
0693       //-------------------------------------
0694       // Region Based Kinematic Cuts: pt of the diphoton system and pt of each photon
0695       //-------------------------------------
0696       if (detector == EcalBarrel) {  // BARREL
0697         //EB region 1
0698         if (etapair <= region1_Barrel_) {  //EB region 1
0699           if (ptmin < selePtGammaBarrel_region1_ || pt_pair < selePtPairBarrel_region1_)
0700             continue;
0701         }
0702         //EB region 2
0703         else {
0704           if (ptmin < selePtGammaBarrel_region2_ || pt_pair < selePtPairBarrel_region2_)
0705             continue;
0706         }
0707       } else {  // ENDCAP
0708         //EE region 1
0709         if (etapair <= region1_EndCap_) {
0710           if (ptmin < selePtGammaEndCap_region1_ || pt_pair < selePtPairEndCap_region1_)
0711             continue;
0712         }
0713         //EE region 2
0714         else if (etapair <= region2_EndCap_) {
0715           if (ptmin < selePtGammaEndCap_region2_ || pt_pair < selePtPairEndCap_region2_)
0716             continue;
0717         }
0718         //EE region 3
0719         else {
0720           if (ptmin < selePtGammaEndCap_region3_ || pt_pair < selePtPairEndCap_region3_)
0721             continue;
0722           if (pt_pair > selePtPairMaxEndCap_region3_)
0723             continue;  // there is also a possible maximum pt for the very forward region
0724         }
0725       }
0726 
0727       /// Mass window Cut
0728       if (m_pair < m_minCut || m_pair > m_maxCut)
0729         continue;
0730 
0731       //// Loop on cluster to measure isolation:
0732       vector<int> IsoClus;
0733       IsoClus.push_back(ind1);  //first two are good clusters
0734       IsoClus.push_back(ind2);
0735 
0736       float Iso = 0;
0737       for (auto it_bc3 = clusterCollection->begin(); it_bc3 != clusterCollection->end(); it_bc3++) {
0738         if (it_bc3->seed() == it_bc->seed() || it_bc3->seed() == it_bc2->seed())
0739           continue;
0740         float drcl = std::hypot(eta_pair - it_bc3->eta(), deltaPhi(phi_pair, it_bc3->phi()));
0741         float dretacl = fabs(eta_pair - it_bc3->eta());
0742         if (drcl > isoBeltdrCut || dretacl > isoBeltdetaCut)
0743           continue;
0744         float pt3 = it_bc3->energy() * sin(it_bc3->position().theta());
0745         if (pt3 < ptminforIso)
0746           continue;
0747         Iso += pt3;
0748         int ind3 = int(it_bc3 - clusterCollection->begin());  /// remember which Iso cluster used
0749         IsoClus.push_back(ind3);
0750       }
0751 
0752       //-------------------------------------
0753       // Region Based Nxtal Cut: Nxtal of the clusters
0754       //-------------------------------------
0755 
0756       uint32_t Nxtal_pho1 = it_bc->size();
0757       uint32_t Nxtal_pho2 = it_bc2->size();
0758 
0759       // BARREL
0760       if (detector == EcalBarrel) {
0761         //EB region 1
0762         if (etapair <= region1_Barrel_) {  //EB region 1
0763           if (Nxtal_pho1 < seleNxtalBarrel_region1_ || Nxtal_pho2 < seleNxtalBarrel_region1_)
0764             continue;
0765         }
0766         //EB region 2
0767         else {
0768           if (Nxtal_pho1 < seleNxtalBarrel_region2_ || Nxtal_pho2 < seleNxtalBarrel_region2_)
0769             continue;
0770         }
0771       }
0772       // ENDCAP
0773       else {
0774         //EE region 1
0775         if (etapair <= region1_EndCap_) {
0776           if (Nxtal_pho1 < seleNxtalEndCap_region1_ || Nxtal_pho2 < seleNxtalEndCap_region1_)
0777             continue;
0778         }
0779         //EE region 2
0780         else if (etapair <= region2_EndCap_) {
0781           if (Nxtal_pho1 < seleNxtalEndCap_region2_ || Nxtal_pho2 < seleNxtalEndCap_region2_)
0782             continue;
0783         }
0784         //EE region 3
0785         else {
0786           if (Nxtal_pho1 < seleNxtalEndCap_region3_ || Nxtal_pho2 < seleNxtalEndCap_region3_)
0787             continue;
0788         }
0789       }
0790 
0791       //-------------------------------------
0792       // Region Based Isolation Cut: pt of the diphoton system and pt of each photon
0793       //-------------------------------------
0794       float iso_val = Iso / pt_pair;
0795       // BARREL
0796       if (detector == EcalBarrel) {
0797         //EB region 1
0798         if (etapair <= region1_Barrel_) {  //EB region 1
0799           if (iso_val > seleIsoBarrel_region1_)
0800             continue;
0801         }
0802         //EB region 2
0803         else {
0804           if (iso_val > seleIsoBarrel_region2_)
0805             continue;
0806         }
0807       }
0808       // ENDCAP
0809       else {
0810         //EE region 1
0811         if (etapair <= region1_EndCap_) {
0812           if (iso_val > seleIsoEndCap_region1_)
0813             continue;
0814         }
0815         //EE region 2
0816         else if (etapair <= region2_EndCap_) {
0817           if (iso_val > seleIsoEndCap_region2_)
0818             continue;
0819         }
0820         //EE region 3
0821         else {
0822           if (iso_val > seleIsoEndCap_region3_)
0823             continue;
0824         }
0825       }
0826       //-------------------------------------
0827       //if(Iso/pt_pair > isoCut) continue;  //old non-regional cut
0828 
0829       //-------------------------------------
0830       // Region based ShowerShape Cut: S4S9 Cut with possible S9S25 configuration
0831       //-------------------------------------
0832 
0833       bool failShowerShape = false;
0834       for (int n = 0; n < 2; n++) {
0835         int ind = IsoClus[n];
0836         auto it_bc3 = clusterCollection->begin() + ind;
0837         auto itmap = passShowerShape_clus.find(ind);
0838 
0839         if (itmap != passShowerShape_clus.end()) {  // if we havent reached the end
0840           if (itmap->second == false) {
0841             failShowerShape = true;
0842             n = 2;  //exit the loop
0843           }
0844         } else {
0845           vector<EcalRecHit> RecHitsCluster_5x5;
0846           float res[3];
0847 
0848           //build the shower shape variables
0849           calcShowerShape(*it_bc3, channelStatus, hitCollection, topology_p, store5x5, RecHitsCluster_5x5, res);
0850           float s4s9 = res[1] > 0 ? res[0] / res[1] : -1;
0851           float s9s25 = res[2] > 0 ? res[1] / res[2] : -1;
0852 
0853           bool passed = false;
0854 
0855           // BARREL
0856           if (detector == EcalBarrel) {
0857             //EB region 1
0858             if (etapair <= region1_Barrel_) {
0859               passed = s4s9 > seleS4S9GammaBarrel_region1_;
0860             }
0861             //EB region 2
0862             else {
0863               passed = s4s9 > seleS4S9GammaBarrel_region2_;
0864             }
0865           }
0866           // ENDCAP
0867           else {
0868             //EE region 1
0869             if (etapair <= region1_EndCap_) {
0870               passed = s4s9 > seleS4S9GammaEndCap_region1_;
0871             }
0872             //EE region 2
0873             else if (etapair <= region2_EndCap_) {
0874               passed = s4s9 > seleS4S9GammaEndCap_region2_;
0875             }
0876             //EE region 3
0877             else {
0878               passed = s4s9 > seleS4S9GammaEndCap_region3_;
0879             }
0880           }
0881 
0882           //apply the S9S25 cut as well
0883           passed = passed && (s9s25 > s9s25Cut);
0884           passShowerShape_clus.insert(pair<int, bool>(ind, passed));
0885 
0886           if (!passed) {
0887             failShowerShape = true;
0888             n = 2;  //exit the loop
0889           } else {
0890             RecHits5x5_clus.insert(pair<int, vector<EcalRecHit> >(ind, RecHitsCluster_5x5));
0891           }
0892         }
0893       }
0894 
0895       if (failShowerShape == true)
0896         continue;  //if any of the two clusters fail shower shape
0897 
0898       ///Save two good clusters' index( plus Iso cluster )  if not yet saved
0899       for (unsigned int iii = 0; iii < IsoClus.size(); iii++) {
0900         int ind = IsoClus[iii];
0901 
0902         if (iii < 2) {
0903           auto it = find(indCandClus.begin(), indCandClus.end(), ind);  ///good cluster
0904           if (it == indCandClus.end())
0905             indCandClus.push_back(ind);
0906           else
0907             continue;
0908         } else {
0909           auto it = find(indIsoClus.begin(), indIsoClus.end(), ind);  //iso cluster
0910           if (it == indIsoClus.end())
0911             indIsoClus.push_back(ind);
0912           else
0913             continue;
0914         }
0915 
0916         auto it = find(indClusSelected.begin(), indClusSelected.end(), ind);
0917         if (it != indClusSelected.end())
0918           continue;
0919         indClusSelected.push_back(ind);
0920       }
0921     }
0922   }
0923 }
0924 
0925 void HLTRegionalEcalResonanceFilter::convxtalid(Int_t &nphi, Int_t &neta) {
0926   // Barrel only
0927   // Output nphi 0...359; neta 0...84; nside=+1 (for eta>0), or 0 (for eta<0).
0928   // neta will be [-85,-1] , or [0,84], the minus sign indicates the z<0 side.
0929 
0930   if (neta > 0)
0931     neta -= 1;
0932   if (nphi > 359)
0933     nphi = nphi - 360;
0934 
0935   // final check
0936   if (nphi > 359 || nphi < 0 || neta < -85 || neta > 84) {
0937     LogError("") << " unexpected fatal error in HLTEcalResonanceFilter::convxtalid " << nphi << " " << neta << " "
0938                  << std::endl;
0939   }
0940 }  //end of convxtalid
0941 
0942 int HLTRegionalEcalResonanceFilter::diff_neta_s(Int_t neta1, Int_t neta2) {
0943   Int_t mdiff;
0944   mdiff = (neta1 - neta2);
0945   return mdiff;
0946 }
0947 
0948 // Calculate the distance in xtals taking into account the periodicity of the Barrel
0949 int HLTRegionalEcalResonanceFilter::diff_nphi_s(Int_t nphi1, Int_t nphi2) {
0950   Int_t mdiff;
0951   if (std::abs(nphi1 - nphi2) < (360 - std::abs(nphi1 - nphi2))) {
0952     mdiff = nphi1 - nphi2;
0953   } else {
0954     mdiff = 360 - std::abs(nphi1 - nphi2);
0955     if (nphi1 > nphi2)
0956       mdiff = -mdiff;
0957   }
0958   return mdiff;
0959 }
0960 
0961 void HLTRegionalEcalResonanceFilter::calcShowerShape(const reco::BasicCluster &bc,
0962                                                      const EcalChannelStatus &channelStatus,
0963                                                      const EcalRecHitCollection *recHits,
0964                                                      const CaloSubdetectorTopology *topology_p,
0965                                                      bool calc5x5,
0966                                                      vector<EcalRecHit> &rechit5x5,
0967                                                      float res[]) {
0968   const std::vector<std::pair<DetId, float> > &vid = bc.hitsAndFractions();
0969 
0970   float e2x2[4] = {0};
0971   float e3x3 = 0;
0972   float e5x5 = 0;
0973   int seedx;
0974   int seedy;
0975   DetId seedId = bc.seed();
0976 
0977   bool InBarrel = true;
0978   if (seedId.subdetId() == EcalBarrel) {
0979     EBDetId ebd = EBDetId(seedId);
0980     seedx = ebd.ieta();
0981     seedy = ebd.iphi();
0982     convxtalid(seedy, seedx);
0983   } else {
0984     InBarrel = false;
0985     EEDetId eed = EEDetId(seedId);
0986     seedx = eed.ix();
0987     seedy = eed.iy();
0988   }
0989   int x, y, dx, dy;
0990 
0991   if (calc5x5) {
0992     rechit5x5.clear();
0993     std::vector<DetId> clus_v5x5;
0994     clus_v5x5 = topology_p->getWindow(seedId, 5, 5);
0995 
0996     for (std::vector<DetId>::const_iterator idItr = clus_v5x5.begin(); idItr != clus_v5x5.end(); idItr++) {
0997       DetId ed = *idItr;
0998       if (InBarrel == true) {
0999         EBDetId ebd = EBDetId(ed);
1000         x = ebd.ieta();
1001         y = ebd.iphi();
1002         convxtalid(y, x);
1003         dx = diff_neta_s(x, seedx);
1004         dy = diff_nphi_s(y, seedy);
1005       } else {
1006         EEDetId eed = EEDetId(ed);
1007         x = eed.ix();
1008         y = eed.iy();
1009         dx = x - seedx;
1010         dy = y - seedy;
1011       }
1012       auto rit = recHits->find(ed);
1013       if (rit == recHits->end())
1014         continue;
1015       if (!checkStatusOfEcalRecHit(channelStatus, *rit))
1016         continue;
1017 
1018       float energy = (*rit).energy();
1019       e5x5 += energy;
1020 
1021       auto idItrF = std::find(vid.begin(), vid.end(), std::make_pair(ed, 1.0f));  ///has to add "f", make it float
1022       if (idItrF == vid.end()) {  ///only store those not belonging to this cluster
1023         rechit5x5.push_back(*rit);
1024       } else {  /// S4, S9 are defined inside the cluster, the same as below.
1025         if (std::abs(dx) <= 1 && std::abs(dy) <= 1) {
1026           if (dx <= 0 && dy <= 0)
1027             e2x2[0] += energy;
1028           if (dx >= 0 && dy <= 0)
1029             e2x2[1] += energy;
1030           if (dx <= 0 && dy >= 0)
1031             e2x2[2] += energy;
1032           if (dx >= 0 && dy >= 0)
1033             e2x2[3] += energy;
1034           e3x3 += energy;
1035         }
1036       }
1037     }
1038 
1039   } else {
1040     for (auto const &idItr : vid) {
1041       DetId ed = idItr.first;
1042       if (InBarrel == true) {
1043         EBDetId ebd = EBDetId(ed);
1044         x = ebd.ieta();
1045         y = ebd.iphi();
1046         convxtalid(y, x);
1047         dx = diff_neta_s(x, seedx);
1048         dy = diff_nphi_s(y, seedy);
1049       } else {
1050         EEDetId eed = EEDetId(ed);
1051         x = eed.ix();
1052         y = eed.iy();
1053         dx = x - seedx;
1054         dy = y - seedy;
1055       }
1056       auto rit = recHits->find(ed);
1057       if (rit == recHits->end()) {
1058         continue;
1059       }
1060 
1061       float energy = (*rit).energy();
1062       if (std::abs(dx) <= 1 && std::abs(dy) <= 1) {
1063         if (dx <= 0 && dy <= 0)
1064           e2x2[0] += energy;
1065         if (dx >= 0 && dy <= 0)
1066           e2x2[1] += energy;
1067         if (dx <= 0 && dy >= 0)
1068           e2x2[2] += energy;
1069         if (dx >= 0 && dy >= 0)
1070           e2x2[3] += energy;
1071         e3x3 += energy;
1072       }
1073     }
1074     e5x5 = e3x3;  ///if not asked to calculte 5x5, then just make e5x5 = e3x3
1075   }
1076 
1077   ///e2x2
1078   res[0] = *max_element(e2x2, e2x2 + 4);
1079   res[1] = e3x3;
1080   res[2] = e5x5;
1081 }
1082 
1083 void HLTRegionalEcalResonanceFilter::calcPaircluster(const reco::BasicCluster &bc1,
1084                                                      const reco::BasicCluster &bc2,
1085                                                      float &m_pair,
1086                                                      float &pt_pair,
1087                                                      float &eta_pair,
1088                                                      float &phi_pair) {
1089   // use XYZVector instead of TLorentzVector to make things faster (and initialize with cartesian coordinates).
1090   // We are interested in the momentum vector:  however, we start from cartesian coordinates to get the vector direction,
1091   // then we set the vector's magnitude to obtain momentum coordinates. The magnitude we set is equal to the particle's energy.
1092   // We can do this because, assuming massless particles (or negligible mass), the magnitude of the momentum vector is given by the energy.
1093 
1094   math::XYZVector v1(bc1.position());
1095   float en1 = bc1.energy();
1096   float scaleFactor = en1 / v1.R();  // XYZVector::R() returns sqrt(Mag2()), where Mag2()= fx*fx + fy*fy + fz*fz
1097   // 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)
1098   v1 *= scaleFactor;
1099 
1100   // 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.
1101   math::XYZVector vsum(bc2.position());
1102   // define energy sum initializing it to energy2, so that we can use it before summing energy1
1103   float energysum = bc2.energy();
1104   scaleFactor = energysum / vsum.R();
1105   vsum *= scaleFactor;
1106 
1107   vsum += v1;
1108   // now sum the energy of the second basic cluster to get total energy
1109   energysum += en1;
1110 
1111   // finally, assign values
1112   m_pair = sqrt(energysum * energysum - vsum.Mag2());  // M_pi0 = sqrt(E_pi0^2 - |p_pi0|^2)
1113   pt_pair = vsum.Rho();  // Rho method is the equivalent of Pt: returns sqrt( fx*fx + fy*fy )
1114   eta_pair = vsum.Eta();
1115   phi_pair = vsum.Phi();
1116 }
1117 
1118 bool HLTRegionalEcalResonanceFilter::checkStatusOfEcalRecHit(const EcalChannelStatus &channelStatus,
1119                                                              const EcalRecHit &rh) {
1120   if (useRecoFlag_) {  ///from recoFlag()
1121     int flag = rh.recoFlag();
1122     if (flagLevelRecHitsToUse_ == 0) {  ///good
1123       if (flag != 0)
1124         return false;
1125     } else if (flagLevelRecHitsToUse_ == 1) {  ///good || PoorCalib
1126       if (flag != 0 && flag != 4)
1127         return false;
1128     } else if (flagLevelRecHitsToUse_ == 2) {  ///good || PoorCalib || LeadingEdgeRecovered || kNeighboursRecovered,
1129       if (flag != 0 && flag != 4 && flag != 6 && flag != 7)
1130         return false;
1131     }
1132   }
1133   if (useDBStatus_) {  //// from DB
1134     int status = int(channelStatus[rh.id().rawId()].getStatusCode());
1135     if (status > statusLevelRecHitsToUse_)
1136       return false;
1137   }
1138 
1139   return true;
1140 }
1141 
1142 void HLTRegionalEcalResonanceFilter::makeClusterES(
1143     float x, float y, float z, const CaloSubdetectorGeometry *geometry_es, const CaloSubdetectorTopology *topology_es) {
1144   ///get assosicated ES clusters of this endcap cluster
1145   const GlobalPoint point(x, y, z);
1146   DetId tmp1 = (dynamic_cast<const EcalPreshowerGeometry *>(geometry_es))->getClosestCellInPlane(point, 1);
1147   DetId tmp2 = (dynamic_cast<const EcalPreshowerGeometry *>(geometry_es))->getClosestCellInPlane(point, 2);
1148   ESDetId strip1 = (tmp1 == DetId(0)) ? ESDetId(0) : ESDetId(tmp1);
1149   ESDetId strip2 = (tmp2 == DetId(0)) ? ESDetId(0) : ESDetId(tmp2);
1150 
1151   // Get ES clusters (found by the PreshSeeded algorithm) associated with a given EE cluster.
1152   for (int i2 = 0; i2 < preshNclust_; i2++) {
1153     reco::PreshowerCluster cl1 =
1154         presh_algo_->makeOneCluster(strip1, &m_used_strips, &m_esrechit_map, geometry_es, topology_es);
1155     reco::PreshowerCluster cl2 =
1156         presh_algo_->makeOneCluster(strip2, &m_used_strips, &m_esrechit_map, geometry_es, topology_es);
1157   }  // end of cycle over ES clusters
1158 }
1159 
1160 // declare this class as a framework plugin
1161 #include "FWCore/Framework/interface/MakerMacros.h"
1162 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!