EgammaClusterProducers/src/UncleanSCRecoveryProducer.cc

0001 /*
0002 UncleanSCRecoveryProducer:
0003 ^^^^^^^^^^^^^^^^^^^^^^^^^^
0004
0005 takes the collections of flagged clean and unclean only SC
0006 (this is the output of UnifiedSCCollectionProducer) and
0007 recovers the original collection of unclean SC.
0008
0009 18 Aug 2010
0010 Nikolaos Rompotis and Chris Seez  - Imperial College London
0011 many thanks to David Wardrope, Shahram Rahatlou and Federico Ferri
0012 */
0013
0014 #include "DataFormats/CaloRecHit/interface/CaloCluster.h"
0015 #include "DataFormats/Common/interface/Handle.h"
0016 #include "DataFormats/EcalDetId/interface/EBDetId.h"
0017 #include "DataFormats/EcalRecHit/interface/EcalRecHit.h"
0018 #include "DataFormats/EcalRecHit/interface/EcalRecHitCollections.h"
0019 #include "DataFormats/EgammaReco/interface/BasicCluster.h"
0020 #include "DataFormats/EgammaReco/interface/BasicClusterFwd.h"
0021 #include "DataFormats/EgammaReco/interface/SuperCluster.h"
0022 #include "DataFormats/EgammaReco/interface/SuperClusterFwd.h"
0023 #include "FWCore/Framework/interface/ESHandle.h"
0024 #include "FWCore/Framework/interface/Event.h"
0025 #include "FWCore/Framework/interface/EventSetup.h"
0026 #include "FWCore/Framework/interface/global/EDProducer.h"
0027 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0028 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0029
0030 #include <iostream>
0031 #include <memory>
0032 #include <vector>
0033
0034 class UncleanSCRecoveryProducer : public edm::global::EDProducer<> {
0035 public:
0036   UncleanSCRecoveryProducer(const edm::ParameterSet& ps);
0037
0038   void produce(edm::StreamID, edm::Event&, const edm::EventSetup&) const override;
0039
0040 private:
0041   // the clean collection
0042   const edm::EDGetTokenT<reco::BasicClusterCollection> cleanBcCollection_;
0043   const edm::EDGetTokenT<reco::SuperClusterCollection> cleanScCollection_;
0044   // the uncleaned collection
0045   const edm::EDGetTokenT<reco::BasicClusterCollection> uncleanBcCollection_;
0046   const edm::EDGetTokenT<reco::SuperClusterCollection> uncleanScCollection_;
0047   // the names of the products to be produced:
0048   const std::string bcCollection_;
0049   const std::string scCollection_;
0050 };
0051
0052 #include "FWCore/Framework/interface/MakerMacros.h"
0053 DEFINE_FWK_MODULE(UncleanSCRecoveryProducer);
0054
0055 UncleanSCRecoveryProducer::UncleanSCRecoveryProducer(const edm::ParameterSet& ps)
0056     : cleanBcCollection_(consumes<reco::BasicClusterCollection>(ps.getParameter<edm::InputTag>("cleanBcCollection"))),
0057       cleanScCollection_(consumes<reco::SuperClusterCollection>(ps.getParameter<edm::InputTag>("cleanScCollection"))),
0058       uncleanBcCollection_(
0059           consumes<reco::BasicClusterCollection>(ps.getParameter<edm::InputTag>("uncleanBcCollection"))),
0060       uncleanScCollection_(
0061           consumes<reco::SuperClusterCollection>(ps.getParameter<edm::InputTag>("uncleanScCollection"))),
0062       bcCollection_(ps.getParameter<std::string>("bcCollection")),
0063       scCollection_(ps.getParameter<std::string>("scCollection")) {
0064   // the products:
0065   produces<reco::BasicClusterCollection>(bcCollection_);
0066   produces<reco::SuperClusterCollection>(scCollection_);
0067 }
0068
0069 void UncleanSCRecoveryProducer::produce(edm::StreamID, edm::Event& evt, const edm::EventSetup& es) const {
0070   // __________________________________________________________________________
0071   //
0072   // cluster collections:
0073   edm::Handle<reco::BasicClusterCollection> pCleanBC;
0074   edm::Handle<reco::SuperClusterCollection> pCleanSC;
0075   //
0076   edm::Handle<reco::BasicClusterCollection> pUncleanBC;
0077   edm::Handle<reco::SuperClusterCollection> pUncleanSC;
0078   // clean collections ________________________________________________________
0079   evt.getByToken(cleanScCollection_, pCleanSC);
0080   const reco::SuperClusterCollection cleanSC = *(pCleanSC.product());
0081
0082   // unclean collections ______________________________________________________
0083   evt.getByToken(uncleanBcCollection_, pUncleanBC);
0084   const reco::BasicClusterCollection uncleanBC = *(pUncleanBC.product());
0085   //
0086   evt.getByToken(uncleanScCollection_, pUncleanSC);
0087   const reco::SuperClusterCollection uncleanSC = *(pUncleanSC.product());
0088   int uncleanSize = pUncleanSC->size();
0089   int cleanSize = pCleanSC->size();
0090
0091   LogTrace("EcalCleaning") << "Size of Clean Collection: " << cleanSize << ", uncleanSize: " << uncleanSize;
0092
0093   // collections are all taken now ____________________________________________
0094   //
0095   // the collections to be produced ___________________________________________
0096   reco::BasicClusterCollection basicClusters;
0097   reco::SuperClusterCollection superClusters;
0098   //
0099   //
0100   // collect all the basic clusters of the SC that belong to the unclean
0101   // collection and put them into the basicClusters vector
0102   // keep the information of which SC they belong to
0103   //
0104   // loop over the unclean sc: all SC will join the new collection
0105   std::vector<std::pair<int, int> > basicClusterOwner;  // counting all
0106
0107   std::vector<DetId> scUncleanSeedDetId;  // counting the unclean
0108   for (int isc = 0; isc < uncleanSize; ++isc) {
0109     const reco::SuperCluster& unsc = uncleanSC[isc];
0110     scUncleanSeedDetId.push_back(unsc.seed()->seed());
0111     reco::CaloCluster_iterator bciter = unsc.clustersBegin();
0112     for (; bciter != unsc.clustersEnd(); ++bciter) {
0113       // the basic clusters
0114       basicClusters.push_back(**bciter);
0115       // index of the unclean SC
0116       basicClusterOwner.push_back(std::make_pair(isc, 0));
0117     }
0118   }
0119   // loop over the clean: only the ones which are in common with the unclean
0120   // are taken into account
0121
0122   std::vector<DetId> scCleanSeedDetId;  // counting the clean
0123   std::vector<int> isToBeKept;
0124   for (int isc = 0; isc < cleanSize; ++isc) {
0125     reco::SuperClusterRef cscRef(pCleanSC, isc);
0126     scCleanSeedDetId.push_back(cscRef->seed()->seed());
0127     for (reco::CaloCluster_iterator bciter = cscRef->clustersBegin(); bciter != cscRef->clustersEnd(); ++bciter) {
0128       // the basic clusters
0129       basicClusters.push_back(**bciter);
0130       // index of the clean SC
0131       basicClusterOwner.push_back(std::make_pair(isc, 1));
0132     }
0133     if (cscRef->isInUnclean())
0134       isToBeKept.push_back(1);
0135     else
0136       isToBeKept.push_back(0);
0137   }
0138   //
0139   // now export the basic clusters into the event and get them back
0140   auto basicClusters_p = std::make_unique<reco::BasicClusterCollection>();
0141   basicClusters_p->assign(basicClusters.begin(), basicClusters.end());
0142   edm::OrphanHandle<reco::BasicClusterCollection> bccHandle = evt.put(std::move(basicClusters_p), bcCollection_);
0143   if (!(bccHandle.isValid())) {
0144     edm::LogWarning("MissingInput") << "could not handle the new BasicClusters!";
0145     return;
0146   }
0147
0148   LogTrace("EcalCleaning") << "Got the BasicClusters from the event again";
0149   int bcSize = bccHandle->size();
0150   //
0151   // now we can create the SC collection
0152   //
0153   // run over the unclean SC: all to be kept here
0154   for (int isc = 0; isc < uncleanSize; ++isc) {
0155     reco::CaloClusterPtrVector clusterPtrVector;
0156     // the seed is the basic cluster with the highest energy
0157     reco::CaloClusterPtr seed;
0158     for (int jbc = 0; jbc < bcSize; ++jbc) {
0159       std::pair<int, int> theBcOwner = basicClusterOwner[jbc];
0160       if (theBcOwner.first == isc && theBcOwner.second == 0) {
0161         reco::CaloClusterPtr currentClu = reco::CaloClusterPtr(bccHandle, jbc);
0162         clusterPtrVector.push_back(currentClu);
0163         if (scUncleanSeedDetId[isc] == currentClu->seed()) {
0164           seed = currentClu;
0165         }
0166       }
0167     }
0168     const reco::SuperCluster& unsc = uncleanSC[isc];
0169     reco::SuperCluster newSC(unsc.energy(), unsc.position(), seed, clusterPtrVector);
0170     newSC.setFlags(reco::CaloCluster::uncleanOnly);
0171     superClusters.push_back(newSC);
0172   }
0173   // run over the clean SC: only those who are in common between the
0174   // clean and unclean collection are kept
0175   for (int isc = 0; isc < cleanSize; ++isc) {
0176     reco::SuperClusterRef cscRef(pCleanSC, isc);
0177     if (not cscRef->isInUnclean())
0178       continue;
0179     reco::CaloClusterPtrVector clusterPtrVector;
0180     // the seed is the basic cluster with the highest energy
0181     reco::CaloClusterPtr seed;
0182     for (int jbc = 0; jbc < bcSize; ++jbc) {
0183       std::pair<int, int> theBcOwner = basicClusterOwner[jbc];
0184       if (theBcOwner.first == isc && theBcOwner.second == 1) {
0185         reco::CaloClusterPtr currentClu = reco::CaloClusterPtr(bccHandle, jbc);
0186         clusterPtrVector.push_back(currentClu);
0187         if (scCleanSeedDetId[isc] == currentClu->seed()) {
0188           seed = currentClu;
0189         }
0190       }
0191     }
0192     reco::SuperCluster newSC(cscRef->energy(), cscRef->position(), seed, clusterPtrVector);
0193     newSC.setFlags(reco::CaloCluster::common);
0194     superClusters.push_back(newSC);
0195   }
0196
0197   auto superClusters_p = std::make_unique<reco::SuperClusterCollection>();
0198   superClusters_p->assign(superClusters.begin(), superClusters.end());
0199
0200   evt.put(std::move(superClusters_p), scCollection_);
0201
0202   LogTrace("EcalCleaning") << "Clusters (Basic/Super) added to the Event! :-)";
0203
0204   // ----- debugging ----------
0205   // print the new collection SC quantities
0206   // print out the clean collection SC
0207   LogTrace("EcalCleaning") << "Clean Collection SC ";
0208   for (int i = 0; i < cleanSize; ++i) {
0209     const reco::SuperCluster& csc = cleanSC[i];
0210     LogTrace("EcalCleaning") << " >>> clean    #" << i << "; Energy: " << csc.energy() << " eta: " << csc.eta()
0211                              << " sc seed detid: " << csc.seed()->seed().rawId();
0212   }
0213   // the unclean SC
0214   LogTrace("EcalCleaning") << "Unclean Collection SC ";
0215   for (int i = 0; i < uncleanSize; ++i) {
0216     const reco::SuperCluster& usc = uncleanSC[i];
0217     LogTrace("EcalCleaning") << " >>> unclean  #" << i << "; Energy: " << usc.energy() << " eta: " << usc.eta()
0218                              << " sc seed detid: " << usc.seed()->seed().rawId();
0219   }
0220   // the new collection
0221   LogTrace("EcalCleaning") << "The new SC clean collection with size " << superClusters.size();
0222   for (unsigned int i = 0; i < superClusters.size(); ++i) {
0223     const reco::SuperCluster nsc = superClusters[i];
0224     LogTrace("EcalCleaning") << " >>> newSC    #" << i << "; Energy: " << nsc.energy() << " eta: " << nsc.eta()
0225                              << " isClean=" << nsc.isInClean() << " isUnclean=" << nsc.isInUnclean()
0226                              << " sc seed detid: " << nsc.seed()->seed().rawId();
0227   }
0228 }