EcalAlCaRecoProducers/plugins/AlCaECALRecHitReducer.cc

0001 #include "Calibration/EcalAlCaRecoProducers/plugins/AlCaECALRecHitReducer.h"
0002 #include "DataFormats/EgammaCandidates/interface/Photon.h"
0003 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0004
0005 #include "DataFormats/EgammaReco/interface/BasicCluster.h"
0006 #include "DataFormats/EcalDetId/interface/EcalSubdetector.h"
0007 #include "FWCore/Framework/interface/ESHandle.h"
0008 #include "FWCore/Framework/interface/MakerMacros.h"
0009
0010 //#define ALLrecHits
0011 //#define DEBUG
0012
0013 //#define QUICK -> if commented loop over the recHits of the SC and add them to the list of recHits to be saved
0014 //                 comment it if you want a faster module but be sure the window is large enough
0015
0016 /** This module reduces the recHitCollections and the caloCalusterCollections in input
0017  * keeping only those associated to the given electrons/photons
0018  */
0019
0020 /// \todo make sure that the new caloClusterCollection has no duplicates
0021
0022 AlCaECALRecHitReducer::AlCaECALRecHitReducer(const edm::ParameterSet& iConfig) {
0023   ebRecHitsToken_ = consumes<EcalRecHitCollection>(iConfig.getParameter<edm::InputTag>("ebRecHitsLabel"));
0024   eeRecHitsToken_ = consumes<EcalRecHitCollection>(iConfig.getParameter<edm::InputTag>("eeRecHitsLabel"));
0025   //  esRecHitsLabel_ = iConfig.getParameter< edm::InputTag > ("esRecHitsLabel");
0026
0027   std::vector<edm::InputTag> srcLabels = iConfig.getParameter<std::vector<edm::InputTag> >("srcLabels");
0028   for (auto inputTag = srcLabels.begin(); inputTag != srcLabels.end(); ++inputTag) {
0029     eleViewTokens_.push_back(consumes<edm::View<reco::RecoCandidate> >(*inputTag));
0030   }
0031
0032   //eleViewToken_ = consumes<edm::View <reco::RecoCandidate> > (iConfig.getParameter< edm::InputTag > ("electronLabel"));
0033   photonToken_ = consumes<reco::PhotonCollection>(iConfig.getParameter<edm::InputTag>("photonLabel"));
0034   EESuperClusterToken_ =
0035       consumes<reco::SuperClusterCollection>(iConfig.getParameter<edm::InputTag>("EESuperClusterCollection"));
0036
0037   caloTopologyToken_ = esConsumes<CaloTopology, CaloTopologyRecord>();
0038
0039   minEta_highEtaSC_ = iConfig.getParameter<double>("minEta_highEtaSC");
0040
0041   alcaBarrelHitsCollection_ = iConfig.getParameter<std::string>("alcaBarrelHitCollection");
0042   alcaEndcapHitsCollection_ = iConfig.getParameter<std::string>("alcaEndcapHitCollection");
0043   alcaCaloClusterCollection_ = iConfig.getParameter<std::string>("alcaCaloClusterCollection");
0044
0045   //  alcaPreshowerHitsCollection_ = iConfig.getParameter<std::string>("alcaPreshowerHitCollection");
0046
0047   etaSize_ = iConfig.getParameter<int>("etaSize");
0048   phiSize_ = iConfig.getParameter<int>("phiSize");
0049   // FIXME: minimum size of etaSize_ and phiSize_
0050   if (phiSize_ % 2 == 0 || etaSize_ % 2 == 0)
0051     edm::LogError("AlCaECALRecHitReducerError") << "Size of eta/phi should be odd numbers";
0052
0053   //  esNstrips_  = iConfig.getParameter<int> ("esNstrips");
0054   //  esNcolumns_ = iConfig.getParameter<int> ("esNcolumns");
0055
0056   //register your products
0057   produces<EBRecHitCollection>(alcaBarrelHitsCollection_);
0058   produces<EERecHitCollection>(alcaEndcapHitsCollection_);
0059   produces<reco::CaloClusterCollection>(alcaCaloClusterCollection_);
0060   //  produces< ESRecHitCollection > (alcaPreshowerHitsCollection_) ;
0061 }
0062
0063 AlCaECALRecHitReducer::~AlCaECALRecHitReducer() {}
0064
0065 // ------------ method called to produce the data  ------------
0066 void AlCaECALRecHitReducer::produce(edm::StreamID, edm::Event& iEvent, const edm::EventSetup& iSetup) const {
0067   using namespace edm;
0068   //using namespace std;
0069
0070   const CaloTopology* caloTopology = &iSetup.getData(caloTopologyToken_);
0071
0072   // Get Photons
0073   Handle<reco::PhotonCollection> phoHandle;
0074   iEvent.getByToken(photonToken_, phoHandle);
0075
0076   // get RecHits
0077   Handle<EBRecHitCollection> barrelRecHitsHandle;
0078   iEvent.getByToken(ebRecHitsToken_, barrelRecHitsHandle);
0079   const EBRecHitCollection* barrelHitsCollection = barrelRecHitsHandle.product();
0080
0081   // get RecHits
0082   Handle<EERecHitCollection> endcapRecHitsHandle;
0083   iEvent.getByToken(eeRecHitsToken_, endcapRecHitsHandle);
0084   const EERecHitCollection* endcapHitsCollection = endcapRecHitsHandle.product();
0085
0086   //   // get ES RecHits
0087   //   Handle<ESRecHitCollection> preshowerRecHitsHandle;
0088   //   iEvent.getByToken(esRecHitsToken_,preshowerRecHitsHandle);
0089
0090   //   const ESRecHitCollection * preshowerHitsCollection = 0 ;
0091   //   if (preshowerIsFull)
0092   //     preshowerHitsCollection = preshowerRecHitsHandle.product () ;
0093
0094   //   // make a vector to store the used ES rechits:
0095   //   set<ESDetId> used_strips;
0096   //   used_strips.clear();
0097
0098   // for Z->ele+SC
0099   Handle<reco::SuperClusterCollection> EESCHandle;
0100   iEvent.getByToken(EESuperClusterToken_, EESCHandle);
0101
0102   //Create empty output collections
0103   auto miniEBRecHitCollection = std::make_unique<EBRecHitCollection>();
0104   auto miniEERecHitCollection = std::make_unique<EERecHitCollection>();
0105   //  auto miniESRecHitCollection = std::make_unique<ESRecHitCollection>();
0106
0107   std::set<DetId> reducedRecHit_EBmap;
0108   std::set<DetId> reducedRecHit_EEmap;
0109
0110   //  std::set< edm::Ref<reco::CaloCluster> > reducedCaloClusters_map;
0111
0112   auto reducedCaloClusterCollection = std::make_unique<reco::CaloClusterCollection>();
0113
0114   //Photons:
0115 #ifdef shervin
0116   for (reco::PhotonCollection::const_iterator phoIt = phoHandle->begin(); phoIt != phoHandle->end(); phoIt++) {
0117     const reco::SuperCluster& sc = *(phoIt->superCluster());
0118
0119     if (phoIt->isEB()) {
0120       AddMiniRecHitCollection(sc, reducedRecHit_EBmap, caloTopology);
0121     } else {  // endcap
0122       AddMiniRecHitCollection(sc, reducedRecHit_EEmap, caloTopology);
0123     }  // end of endcap
0124
0125     /// \todo check if this works when you ask sc->seed(), I suspect that the references have to be updated
0126     reco::CaloCluster_iterator it = sc.clustersBegin();
0127     reco::CaloCluster_iterator itend = sc.clustersEnd();
0128     for (; it != itend; ++it) {
0129       reco::CaloCluster caloClus(**it);
0130       reducedCaloClusterCollection->push_back(caloClus);
0131     }
0132   }
0133 #endif
0134
0135   Handle<edm::View<reco::RecoCandidate> > eleViewHandle;
0136   for (auto iToken = eleViewTokens_.begin(); iToken != eleViewTokens_.end(); iToken++) {
0137     iEvent.getByToken(*iToken, eleViewHandle);
0138
0139     //Electrons:
0140     for (auto eleIt = eleViewHandle->begin(); eleIt != eleViewHandle->end(); eleIt++) {
0141       const reco::SuperCluster& sc = *(eleIt->superCluster());
0142
0143       if (fabs(sc.eta()) < 1.479) {
0144         AddMiniRecHitCollection(sc, reducedRecHit_EBmap, caloTopology);
0145       } else {  // endcap
0146         AddMiniRecHitCollection(sc, reducedRecHit_EEmap, caloTopology);
0147       }  // end of endcap
0148
0149       reco::CaloCluster_iterator it = sc.clustersBegin();
0150       reco::CaloCluster_iterator itend = sc.clustersEnd();
0151       for (; it != itend; ++it) {
0152         reco::CaloCluster caloClus(**it);
0153         reducedCaloClusterCollection->push_back(caloClus);
0154       }
0155     }
0156   }
0157
0158   //saving recHits for highEta SCs for highEta studies
0159   for (reco::SuperClusterCollection::const_iterator SC_iter = EESCHandle->begin(); SC_iter != EESCHandle->end();
0160        SC_iter++) {
0161     if (fabs(SC_iter->eta()) < minEta_highEtaSC_)
0162       continue;
0163     AddMiniRecHitCollection(*SC_iter, reducedRecHit_EEmap, caloTopology);
0164
0165     const reco::SuperCluster& sc = *(SC_iter);
0166     reco::CaloCluster_iterator it = sc.clustersBegin();
0167     reco::CaloCluster_iterator itend = sc.clustersEnd();
0168     for (; it != itend; ++it) {
0169       reco::CaloCluster caloClus(**it);
0170       reducedCaloClusterCollection->push_back(caloClus);
0171     }
0172   }
0173
0174   //------------------------------ fill the alcareco reduced recHit collection
0175   for (std::set<DetId>::const_iterator itr = reducedRecHit_EBmap.begin(); itr != reducedRecHit_EBmap.end(); itr++) {
0176     if (barrelHitsCollection->find(*itr) != barrelHitsCollection->end())
0177       miniEBRecHitCollection->push_back(*(barrelHitsCollection->find(*itr)));
0178   }
0179
0180   for (std::set<DetId>::const_iterator itr = reducedRecHit_EEmap.begin(); itr != reducedRecHit_EEmap.end(); itr++) {
0181     if (endcapHitsCollection->find(*itr) != endcapHitsCollection->end())
0182       miniEERecHitCollection->push_back(*(endcapHitsCollection->find(*itr)));
0183   }
0184
0185   //--------------------------------------- Put selected information in the event
0186   iEvent.put(std::move(miniEBRecHitCollection), alcaBarrelHitsCollection_);
0187   iEvent.put(std::move(miniEERecHitCollection), alcaEndcapHitsCollection_);
0188   //  iEvent.put(std::move(miniESRecHitCollection),alcaPreshowerHitsCollection_ );
0189   iEvent.put(std::move(reducedCaloClusterCollection), alcaCaloClusterCollection_);
0190 }
0191
0192 void AlCaECALRecHitReducer::AddMiniRecHitCollection(const reco::SuperCluster& sc,
0193                                                     std::set<DetId>& reducedRecHitMap,
0194                                                     const CaloTopology* caloTopology) const {
0195   DetId seed = (sc.seed()->seed());
0196   int phiSize = phiSize_, etaSize = etaSize_;
0197   if (seed.subdetId() != EcalBarrel) {  // if not EB, take a square window
0198     etaSize = std::max(phiSize_, etaSize_);
0199     phiSize = etaSize;
0200   }
0201
0202   std::vector<DetId> recHit_window = caloTopology->getWindow(seed, phiSize, etaSize);
0203   for (unsigned int i = 0; i < recHit_window.size(); i++) {
0204     reducedRecHitMap.insert(recHit_window[i]);
0205   }
0206
0207   const std::vector<std::pair<DetId, float> >& scHits = sc.hitsAndFractions();
0208   for (std::vector<std::pair<DetId, float> >::const_iterator scHit_itr = scHits.begin(); scHit_itr != scHits.end();
0209        scHit_itr++) {
0210     // the map fills just one time (avoiding double insert of recHits)
0211     reducedRecHitMap.insert(scHit_itr->first);
0212   }
0213
0214   return;
0215 }
0216
0217 DEFINE_FWK_MODULE(AlCaECALRecHitReducer);