EgammaClusterProducers/src/CleanAndMergeProducer.cc

0001 #include "CondFormats/DataRecord/interface/EcalChannelStatusRcd.h"
0002 #include "CondFormats/EcalObjects/interface/EcalChannelStatus.h"
0003 #include "DataFormats/CaloRecHit/interface/CaloCluster.h"
0004 #include "DataFormats/Common/interface/Handle.h"
0005 #include "DataFormats/EcalDetId/interface/EBDetId.h"
0006 #include "DataFormats/EcalRecHit/interface/EcalRecHit.h"
0007 #include "DataFormats/EcalRecHit/interface/EcalRecHitCollections.h"
0008 #include "DataFormats/EgammaReco/interface/BasicCluster.h"
0009 #include "DataFormats/EgammaReco/interface/BasicClusterFwd.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 "Geometry/CaloGeometry/interface/CaloCellGeometry.h"
0019 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
0020 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
0021 #include "Geometry/CaloTopology/interface/EcalBarrelTopology.h"
0022 #include "Geometry/CaloTopology/interface/EcalEndcapTopology.h"
0023 #include "Geometry/CaloTopology/interface/EcalPreshowerTopology.h"
0024 #include "Geometry/Records/interface/CaloGeometryRecord.h"
0025 #include "RecoEcal/EgammaCoreTools/interface/ClusterShapeAlgo.h"
0026 #include "RecoEcal/EgammaCoreTools/interface/PositionCalc.h"
0027
0028 #include <iostream>
0029 #include <memory>
0030 #include <vector>
0031
0032 class CleanAndMergeProducer : public edm::stream::EDProducer<> {
0033 public:
0034   CleanAndMergeProducer(const edm::ParameterSet& ps);
0035
0036   ~CleanAndMergeProducer() override;
0037
0038   void produce(edm::Event&, const edm::EventSetup&) override;
0039
0040 private:
0041   edm::EDGetTokenT<reco::SuperClusterCollection> cleanScToken_;
0042   edm::EDGetTokenT<reco::SuperClusterCollection> uncleanScToken_;
0043
0044   // the names of the products to be produced:
0045   std::string bcCollection_;
0046   std::string scCollection_;
0047   std::string refScCollection_;
0048 };
0049
0050 #include "FWCore/Framework/interface/MakerMacros.h"
0051 DEFINE_FWK_MODULE(CleanAndMergeProducer);
0052
0053 /*
0054 CleanAndMergeProducer:
0055 ^^^^^^^^^^^^^^^^^^^^^^
0056
0057 Takes  as  input the cleaned  and the  uncleaned collection of SC
0058 and produces an uncleaned collection of SC without the duplicates
0059 and  a collection of  references  to the cleaned collection of SC
0060
0061 25 June 2010
0062 Nikolaos Rompotis and Chris Seez  - Imperial College London
0063 many thanks to Shahram Rahatlou and Federico Ferri
0064
0065
0066 */
0067
0068 CleanAndMergeProducer::CleanAndMergeProducer(const edm::ParameterSet& ps) {
0069   // get the parameters
0070   // the cleaned collection:
0071   cleanScToken_ = consumes<reco::SuperClusterCollection>(ps.getParameter<edm::InputTag>("cleanScInputTag"));
0072   uncleanScToken_ = consumes<reco::SuperClusterCollection>(ps.getParameter<edm::InputTag>("uncleanScInputTag"));
0073
0074   // the names of the products to be produced:
0075   bcCollection_ = ps.getParameter<std::string>("bcCollection");
0076   scCollection_ = ps.getParameter<std::string>("scCollection");
0077   refScCollection_ = ps.getParameter<std::string>("refScCollection");
0078
0079   std::map<std::string, double> providedParameters;
0080   providedParameters.insert(std::make_pair("LogWeighted", ps.getParameter<bool>("posCalc_logweight")));
0081   providedParameters.insert(std::make_pair("T0_barl", ps.getParameter<double>("posCalc_t0")));
0082   providedParameters.insert(std::make_pair("W0", ps.getParameter<double>("posCalc_w0")));
0083   providedParameters.insert(std::make_pair("X0", ps.getParameter<double>("posCalc_x0")));
0084
0085   // the products:
0086   produces<reco::BasicClusterCollection>(bcCollection_);
0087   produces<reco::SuperClusterCollection>(scCollection_);
0088   produces<reco::SuperClusterRefVector>(refScCollection_);
0089 }
0090
0091 CleanAndMergeProducer::~CleanAndMergeProducer() { ; }
0092
0093 void CleanAndMergeProducer::produce(edm::Event& evt, const edm::EventSetup& es) {
0094   // get the input collections
0095   // ______________________________________________________________________________________
0096   //
0097   // cluster collections:
0098
0099   edm::Handle<reco::SuperClusterCollection> pCleanSC;
0100   edm::Handle<reco::SuperClusterCollection> pUncleanSC;
0101
0102   evt.getByToken(cleanScToken_, pCleanSC);
0103   evt.getByToken(uncleanScToken_, pUncleanSC);
0104
0105   //
0106   // collections are all taken now _____________________________________________________
0107   //
0108   //
0109   // the collections to be produced:
0110   reco::BasicClusterCollection basicClusters;
0111   reco::SuperClusterCollection superClusters;
0112   reco::SuperClusterRefVector* scRefs = new reco::SuperClusterRefVector;
0113
0114   //
0115   // run over the uncleaned SC and check how many of them are matched to the cleaned ones
0116   // if you find a matched one, create a reference to the cleaned collection and store it
0117   // if you find an unmatched one, then keep all its basic clusters in the basic Clusters
0118   // vector
0119   int uncleanSize = pUncleanSC->size();
0120   int cleanSize = pCleanSC->size();
0121
0122   LogTrace("EcalCleaning") << "Size of Clean Collection: " << cleanSize << ", uncleanSize: " << uncleanSize
0123                            << std::endl;
0124   //
0125   // keep whether the SC in unique in the uncleaned collection
0126   std::vector<int> isUncleanOnly;      // 1 if unique in the uncleaned
0127   std::vector<int> basicClusterOwner;  // contains the index of the SC that owns that BS
0128   std::vector<int> isSeed;             // if this basic cluster is a seed it is 1
0129   for (int isc = 0; isc < uncleanSize; ++isc) {
0130     reco::SuperClusterRef unscRef(pUncleanSC, isc);
0131     const std::vector<std::pair<DetId, float> >& uhits = unscRef->hitsAndFractions();
0132     int uhitsSize = uhits.size();
0133     bool foundTheSame = false;
0134     for (int jsc = 0; jsc < cleanSize; ++jsc) {  // loop over the cleaned SC
0135       reco::SuperClusterRef cscRef(pCleanSC, jsc);
0136       const std::vector<std::pair<DetId, float> >& chits = cscRef->hitsAndFractions();
0137       int chitsSize = chits.size();
0138       foundTheSame = true;
0139       if (unscRef->seed()->seed() == cscRef->seed()->seed() && chitsSize == uhitsSize) {
0140         // if the clusters are exactly the same then because the clustering
0141         // algorithm works in a deterministic way, the order of the rechits
0142         // will be the same
0143         for (int i = 0; i < chitsSize; ++i) {
0144           if (uhits[i].first != chits[i].first) {
0145             foundTheSame = false;
0146             break;
0147           }
0148         }
0149         if (foundTheSame) {  // ok you have found it!
0150           // make the reference
0151           //scRefs->push_back( edm::Ref<reco::SuperClusterCollection>(pCleanSC, jsc) );
0152           scRefs->push_back(cscRef);
0153           isUncleanOnly.push_back(0);
0154           break;
0155         }
0156       }
0157     }
0158     if (not foundTheSame) {
0159       // mark it as unique in the uncleaned
0160       isUncleanOnly.push_back(1);
0161       // keep all its basic clusters
0162       for (reco::CaloCluster_iterator bciter = unscRef->clustersBegin(); bciter != unscRef->clustersEnd(); ++bciter) {
0163         // the basic clusters
0164         basicClusters.push_back(**bciter);
0165         basicClusterOwner.push_back(isc);
0166       }
0167     }
0168   }
0169   int bcSize = basicClusters.size();
0170
0171   LogDebug("EcalCleaning") << "Found cleaned SC: " << cleanSize << " uncleaned SC: " << uncleanSize << " from which "
0172                            << scRefs->size() << " will become refs to the cleaned collection";
0173
0174   // now you have the collection of basic clusters of the SC to be remain in the
0175   // in the clean collection, export them to the event
0176   // you will need to reread them later in order to make correctly the refs to the SC
0177   auto basicClusters_p = std::make_unique<reco::BasicClusterCollection>();
0178   basicClusters_p->assign(basicClusters.begin(), basicClusters.end());
0179   edm::OrphanHandle<reco::BasicClusterCollection> bccHandle = evt.put(std::move(basicClusters_p), bcCollection_);
0180   if (!(bccHandle.isValid())) {
0181     edm::LogWarning("MissingInput") << "could not get a handle on the BasicClusterCollection!" << std::endl;
0182     return;
0183   }
0184
0185   LogDebug("EcalCleaning") << "Got the BasicClusters from the event again";
0186   // now you have to create again your superclusters
0187   // you run over the uncleaned SC, but now you know which of them are
0188   // the ones that are needed and which are their basic clusters
0189   for (int isc = 0; isc < uncleanSize; ++isc) {
0190     if (isUncleanOnly[isc] == 1) {  // look for sc that are unique in the unclean collection
0191       // make the basic cluster collection
0192       reco::CaloClusterPtrVector clusterPtrVector;
0193       reco::CaloClusterPtr seed;  // the seed is the basic cluster with the highest energy
0194       double energy = -1;
0195       for (int jbc = 0; jbc < bcSize; ++jbc) {
0196         if (basicClusterOwner[jbc] == isc) {
0197           reco::CaloClusterPtr currentClu = reco::CaloClusterPtr(bccHandle, jbc);
0198           clusterPtrVector.push_back(currentClu);
0199           if (energy < currentClu->energy()) {
0200             energy = currentClu->energy();
0201             seed = currentClu;
0202           }
0203         }
0204       }
0205       reco::SuperClusterRef unscRef(pUncleanSC, isc);
0206       reco::SuperCluster newSC(unscRef->energy(), unscRef->position(), seed, clusterPtrVector);
0207       superClusters.push_back(newSC);
0208     }
0209   }
0210   // export the collection of references to the clean collection
0211   std::unique_ptr<reco::SuperClusterRefVector> scRefs_p(scRefs);
0212   evt.put(std::move(scRefs_p), refScCollection_);
0213
0214   // the collection of basic clusters is already in the event
0215   // the collection of uncleaned SC
0216   auto superClusters_p = std::make_unique<reco::SuperClusterCollection>();
0217   superClusters_p->assign(superClusters.begin(), superClusters.end());
0218   evt.put(std::move(superClusters_p), scCollection_);
0219   LogDebug("EcalCleaning") << "Hybrid Clusters (Basic/Super) added to the Event! :-)";
0220 }