EgammaClusterProducers/src/EcalDigiSelector.cc

0001 #include "DataFormats/DetId/interface/DetId.h"
0002 #include "DataFormats/EcalDetId/interface/EBDetId.h"
0003 #include "DataFormats/EcalDetId/interface/EEDetId.h"
0004 #include "DataFormats/EcalDigi/interface/EBDataFrame.h"
0005 #include "DataFormats/EcalDigi/interface/EEDataFrame.h"
0006 #include "DataFormats/EcalDigi/interface/EcalDataFrame.h"
0007 #include "DataFormats/EcalDigi/interface/EcalDigiCollections.h"
0008 #include "DataFormats/EcalRecHit/interface/EcalRecHitCollections.h"
0009 #include "DataFormats/EgammaReco/interface/BasicCluster.h"
0010 #include "DataFormats/EgammaReco/interface/SuperCluster.h"
0011 #include "DataFormats/EgammaReco/interface/SuperClusterFwd.h"
0012 #include "FWCore/Framework/interface/ESHandle.h"
0013 #include "FWCore/Framework/interface/Event.h"
0014 #include "FWCore/Framework/interface/EventSetup.h"
0015 #include "FWCore/Framework/interface/stream/EDProducer.h"
0016 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0017 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0018 #include "FWCore/Utilities/interface/ESGetToken.h"
0019 #include "Geometry/CaloTopology/interface/CaloTopology.h"
0020 #include "Geometry/Records/interface/CaloTopologyRecord.h"
0021 #include "RecoEcal/EgammaCoreTools/interface/EcalClusterTools.h"
0022
0023 #include <TMath.h>
0024
0025 #include <iostream>
0026 #include <map>
0027 #include <memory>
0028 #include <set>
0029 #include <string>
0030 #include <vector>
0031
0032 class EcalDigiSelector : public edm::stream::EDProducer<> {
0033 public:
0034   EcalDigiSelector(const edm::ParameterSet& ps);
0035
0036   void produce(edm::Event&, const edm::EventSetup&) override;
0037
0038 private:
0039   std::string selectedEcalEBDigiCollection_;
0040   std::string selectedEcalEEDigiCollection_;
0041
0042   edm::EDGetTokenT<reco::SuperClusterCollection> barrelSuperClusterProducer_;
0043   edm::EDGetTokenT<reco::SuperClusterCollection> endcapSuperClusterProducer_;
0044
0045   // input configuration
0046   edm::EDGetTokenT<EcalRecHitCollection> EcalEBRecHitToken_;
0047   edm::EDGetTokenT<EcalRecHitCollection> EcalEERecHitToken_;
0048   edm::EDGetTokenT<EBDigiCollection> EcalEBDigiToken_;
0049   edm::EDGetTokenT<EEDigiCollection> EcalEEDigiToken_;
0050   edm::ESGetToken<CaloTopology, CaloTopologyRecord> caloTopologyToken_;
0051
0052   double cluster_pt_thresh_;
0053   double single_cluster_thresh_;
0054   int nclus_sel_;
0055 };
0056
0057 #include "FWCore/Framework/interface/MakerMacros.h"
0058 DEFINE_FWK_MODULE(EcalDigiSelector);
0059
0060 using namespace reco;
0061 EcalDigiSelector::EcalDigiSelector(const edm::ParameterSet& ps) {
0062   selectedEcalEBDigiCollection_ = ps.getParameter<std::string>("selectedEcalEBDigiCollection");
0063   selectedEcalEEDigiCollection_ = ps.getParameter<std::string>("selectedEcalEEDigiCollection");
0064
0065   barrelSuperClusterProducer_ =
0066       consumes<reco::SuperClusterCollection>(ps.getParameter<edm::InputTag>("barrelSuperClusterProducer"));
0067   endcapSuperClusterProducer_ =
0068       consumes<reco::SuperClusterCollection>(ps.getParameter<edm::InputTag>("endcapSuperClusterProducer"));
0069
0070   EcalEBRecHitToken_ = consumes<EcalRecHitCollection>(ps.getParameter<edm::InputTag>("EcalEBRecHitTag"));
0071   EcalEERecHitToken_ = consumes<EcalRecHitCollection>(ps.getParameter<edm::InputTag>("EcalEERecHitTag"));
0072
0073   EcalEBDigiToken_ = consumes<EBDigiCollection>(ps.getParameter<edm::InputTag>("EcalEBDigiTag"));
0074   EcalEEDigiToken_ = consumes<EEDigiCollection>(ps.getParameter<edm::InputTag>("EcalEEDigiTag"));
0075
0076   caloTopologyToken_ = esConsumes<CaloTopology, CaloTopologyRecord>();
0077
0078   cluster_pt_thresh_ = ps.getParameter<double>("cluster_pt_thresh");
0079   single_cluster_thresh_ = ps.getParameter<double>("single_cluster_thresh");
0080
0081   nclus_sel_ = ps.getParameter<int>("nclus_sel");
0082   produces<EBDigiCollection>(selectedEcalEBDigiCollection_);
0083   produces<EEDigiCollection>(selectedEcalEEDigiCollection_);
0084 }
0085
0086 void EcalDigiSelector::produce(edm::Event& evt, const edm::EventSetup& es) {
0087   //Get BarrelSuperClusters to start.
0088   edm::Handle<reco::SuperClusterCollection> pBarrelSuperClusters;
0089
0090   evt.getByToken(barrelSuperClusterProducer_, pBarrelSuperClusters);
0091
0092   const reco::SuperClusterCollection& BarrelSuperClusters = *pBarrelSuperClusters;
0093   //Got BarrelSuperClusters
0094
0095   //Get BarrelSuperClusters to start.
0096   edm::Handle<reco::SuperClusterCollection> pEndcapSuperClusters;
0097
0098   evt.getByToken(endcapSuperClusterProducer_, pEndcapSuperClusters);
0099
0100   const reco::SuperClusterCollection& EndcapSuperClusters = *pEndcapSuperClusters;
0101   //Got EndcapSuperClusters
0102
0103   reco::SuperClusterCollection saveBarrelSuperClusters;
0104   reco::SuperClusterCollection saveEndcapSuperClusters;
0105   bool meet_single_thresh = false;
0106   //Loop over barrel superclusters, and apply threshold
0107   for (int loop = 0; loop < int(BarrelSuperClusters.size()); loop++) {
0108     const SuperCluster& clus1 = BarrelSuperClusters[loop];
0109     float eta1 = clus1.eta();
0110     float energy1 = clus1.energy();
0111     float theta1 = 2 * atan(exp(-1. * eta1));
0112     float cluspt1 = energy1 * sin(theta1);
0113     if (cluspt1 > cluster_pt_thresh_) {
0114       saveBarrelSuperClusters.push_back(clus1);
0115       if (cluspt1 > single_cluster_thresh_)
0116         meet_single_thresh = true;
0117     }
0118   }
0119
0120   //Loop over endcap superclusters, and apply threshold
0121   for (int loop = 0; loop < int(EndcapSuperClusters.size()); loop++) {
0122     const SuperCluster& clus1 = EndcapSuperClusters[loop];
0123     float eta1 = clus1.eta();
0124     float energy1 = clus1.energy();
0125     float theta1 = 2 * atan(exp(-1. * eta1));
0126     float cluspt1 = energy1 * sin(theta1);
0127     if (cluspt1 > cluster_pt_thresh_) {
0128       saveEndcapSuperClusters.push_back(clus1);
0129       if (cluspt1 > single_cluster_thresh_)
0130         meet_single_thresh = true;
0131     }
0132   }
0133
0134   auto SEBDigiCol = std::make_unique<EBDigiCollection>();
0135   auto SEEDigiCol = std::make_unique<EEDigiCollection>();
0136   int TotClus = saveBarrelSuperClusters.size() + saveEndcapSuperClusters.size();
0137
0138   if (TotClus >= nclus_sel_ || meet_single_thresh) {
0139     if (!saveBarrelSuperClusters.empty()) {
0140       edm::ESHandle<CaloTopology> pTopology = es.getHandle(caloTopologyToken_);
0141       const CaloTopology* topology = pTopology.product();
0142
0143       //get barrel digi collection
0144       edm::Handle<EBDigiCollection> pdigis;
0145       const EBDigiCollection* digis = nullptr;
0146       evt.getByToken(EcalEBDigiToken_, pdigis);
0147       digis = pdigis.product();  // get a ptr to the product
0148
0149       edm::Handle<EcalRecHitCollection> prechits;
0150       const EcalRecHitCollection* rechits = nullptr;
0151       evt.getByToken(EcalEBRecHitToken_, prechits);
0152       rechits = prechits.product();  // get a ptr to the product
0153
0154       if (digis) {
0155         std::vector<DetId> saveTheseDetIds;
0156         //pick out the detids for the 3x3 in each of the selected superclusters
0157         for (int loop = 0; loop < int(saveBarrelSuperClusters.size()); loop++) {
0158           SuperCluster clus1 = saveBarrelSuperClusters[loop];
0159           const CaloClusterPtr& bcref = clus1.seed();
0160           const BasicCluster* bc = bcref.get();
0161           //Get the maximum detid
0162           DetId maxDetId = EcalClusterTools::getMaximum(*bc, rechits).first;
0163           // Loop over the 3x3 array centered on maximum detid
0164           for (DetId detId : CaloRectangleRange(1, maxDetId, *topology))
0165             saveTheseDetIds.push_back(detId);
0166         }
0167         for (int detloop = 0; detloop < int(saveTheseDetIds.size()); ++detloop) {
0168           EBDetId detL = EBDetId(saveTheseDetIds[detloop]);
0169
0170           for (EBDigiCollection::const_iterator blah = digis->begin(); blah != digis->end(); blah++) {
0171             if (detL == blah->id()) {
0172               EBDataFrame myDigi = (*blah);
0173               SEBDigiCol->push_back(detL);
0174
0175               EBDataFrame df(SEBDigiCol->back());
0176               for (int iq = 0; iq < myDigi.size(); ++iq) {
0177                 df.setSample(iq, myDigi.sample(iq).raw());
0178               }
0179               //ebcounter++;
0180             }
0181           }
0182           //if (ebcounter >= int(saveTheseDetIds.size())) break;
0183         }  //loop over dets
0184       }
0185
0186     }  //If barrel superclusters need saving.
0187
0188     if (!saveEndcapSuperClusters.empty()) {
0189       edm::ESHandle<CaloTopology> pTopology = es.getHandle(caloTopologyToken_);
0190       const CaloTopology* topology = pTopology.product();
0191
0192       //Get endcap rec hit collection
0193       //get endcap digi collection
0194       edm::Handle<EEDigiCollection> pdigis;
0195       const EEDigiCollection* digis = nullptr;
0196       evt.getByToken(EcalEEDigiToken_, pdigis);
0197       digis = pdigis.product();  // get a ptr to the product
0198
0199       edm::Handle<EcalRecHitCollection> prechits;
0200       const EcalRecHitCollection* rechits = nullptr;
0201       evt.getByToken(EcalEERecHitToken_, prechits);
0202       rechits = prechits.product();  // get a ptr to the product
0203
0204       if (digis) {
0205         //std::vector<DetId> saveTheseDetIds;
0206         std::set<DetId> saveTheseDetIds;
0207         //pick out the digis for the 3x3 in each of the selected superclusters
0208         for (int loop = 0; loop < int(saveEndcapSuperClusters.size()); loop++) {
0209           SuperCluster clus1 = saveEndcapSuperClusters[loop];
0210           const CaloClusterPtr& bcref = clus1.seed();
0211           const BasicCluster* bc = bcref.get();
0212           //Get the maximum detid
0213           DetId maxDetId = EcalClusterTools::getMaximum(*bc, rechits).first;
0214           // Loop over the 3x3 array centered on maximum detid
0215           for (DetId detId : CaloRectangleRange(1, maxDetId, *topology))
0216             saveTheseDetIds.insert(detId);
0217         }
0218         int eecounter = 0;
0219         for (EEDigiCollection::const_iterator blah = digis->begin(); blah != digis->end(); blah++) {
0220           std::set<DetId>::const_iterator finder = saveTheseDetIds.find(blah->id());
0221           if (finder != saveTheseDetIds.end()) {
0222             EEDetId detL = EEDetId(*finder);
0223
0224             if (detL == blah->id()) {
0225               EEDataFrame myDigi = (*blah);
0226               SEEDigiCol->push_back(detL);
0227               EEDataFrame df(SEEDigiCol->back());
0228               for (int iq = 0; iq < myDigi.size(); ++iq) {
0229                 df.setSample(iq, myDigi.sample(iq).raw());
0230               }
0231               eecounter++;
0232             }
0233           }
0234           if (eecounter >= int(saveTheseDetIds.size()))
0235             break;
0236         }  //loop over digis
0237       }
0238     }  //If endcap superclusters need saving.
0239
0240   }  //If we're actually saving stuff
0241
0242   //Okay, either my collections have been filled with the requisite Digis, or they haven't.
0243
0244   //Empty collection, or full, still put in event.
0245   SEBDigiCol->sort();
0246   SEEDigiCol->sort();
0247   evt.put(std::move(SEBDigiCol), selectedEcalEBDigiCollection_);
0248   evt.put(std::move(SEEDigiCol), selectedEcalEEDigiCollection_);
0249 }