EgammaClusterProducers/src/UnifiedSCCollectionProducer.cc

0001 /*
0002 UnifiedSCCollectionProducer:
0003 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^
0004
0005 Takes  as  input the cleaned  and the  uncleaned collection of SC
0006 and  produces two collections of SC: one with the clean SC, but flagged
0007 such that with the algoID value one can identify the SC that are also
0008 in the unclean collection and a collection with the unclean only SC.
0009 This collection has the algoID enumeration of the SC altered
0010 such that:
0011 flags = 0   (cleanedOnly)     cluster is only in the cleaned collection
0012 flags = 100 (common)          cluster is common in both collections
0013 flags = 200 (uncleanedOnly)   cluster is only in the uncleaned collection
0014
0015 In that way the user can get hold of objects from the
0016 -  cleaned   collection only if they choose flags <  200
0017 -  uncleaned collection only if they choose flags >= 100
0018
0019 18 Aug 2010
0020 Nikolaos Rompotis and Chris Seez  - Imperial College London
0021 many thanks to David Wardrope, Shahram Rahatlou and Federico Ferri
0022 */
0023
0024 #include "CondFormats/DataRecord/interface/EcalChannelStatusRcd.h"
0025 #include "CondFormats/EcalObjects/interface/EcalChannelStatus.h"
0026 #include "DataFormats/CaloRecHit/interface/CaloCluster.h"
0027 #include "DataFormats/Common/interface/Handle.h"
0028 #include "DataFormats/EcalDetId/interface/EBDetId.h"
0029 #include "DataFormats/EcalRecHit/interface/EcalRecHit.h"
0030 #include "DataFormats/EcalRecHit/interface/EcalRecHitCollections.h"
0031 #include "DataFormats/EgammaReco/interface/BasicCluster.h"
0032 #include "DataFormats/EgammaReco/interface/BasicClusterFwd.h"
0033 #include "DataFormats/EgammaReco/interface/SuperCluster.h"
0034 #include "DataFormats/EgammaReco/interface/SuperClusterFwd.h"
0035 #include "FWCore/Framework/interface/ESHandle.h"
0036 #include "FWCore/Framework/interface/Event.h"
0037 #include "FWCore/Framework/interface/EventSetup.h"
0038 #include "FWCore/Framework/interface/stream/EDProducer.h"
0039 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0040 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0041 #include "FWCore/Utilities/interface/InputTag.h"
0042 #include "RecoEcal/EgammaCoreTools/interface/PositionCalc.h"
0043
0044 #include <iostream>
0045 #include <memory>
0046 #include <vector>
0047
0048 class UnifiedSCCollectionProducer : public edm::stream::EDProducer<> {
0049 public:
0050   UnifiedSCCollectionProducer(const edm::ParameterSet& ps);
0051
0052   void produce(edm::Event&, const edm::EventSetup&) override;
0053
0054 private:
0055   // the clean collection
0056   edm::EDGetTokenT<reco::BasicClusterCollection> cleanBcCollection_;
0057   edm::EDGetTokenT<reco::SuperClusterCollection> cleanScCollection_;
0058   // the uncleaned collection
0059   edm::EDGetTokenT<reco::BasicClusterCollection> uncleanBcCollection_;
0060   edm::EDGetTokenT<reco::SuperClusterCollection> uncleanScCollection_;
0061
0062   // the names of the products to be produced:
0063   std::string bcCollection_;
0064   std::string scCollection_;
0065   std::string bcCollectionUncleanOnly_;
0066   std::string scCollectionUncleanOnly_;
0067 };
0068
0069 #include "FWCore/Framework/interface/MakerMacros.h"
0070 DEFINE_FWK_MODULE(UnifiedSCCollectionProducer);
0071
0072 UnifiedSCCollectionProducer::UnifiedSCCollectionProducer(const edm::ParameterSet& ps) {
0073   using reco::BasicClusterCollection;
0074   using reco::SuperClusterCollection;
0075   // get the parameters
0076   // the cleaned collection:
0077   cleanBcCollection_ = consumes<BasicClusterCollection>(ps.getParameter<edm::InputTag>("cleanBcCollection"));
0078   cleanScCollection_ = consumes<SuperClusterCollection>(ps.getParameter<edm::InputTag>("cleanScCollection"));
0079
0080   // the uncleaned collection
0081   uncleanBcCollection_ = consumes<BasicClusterCollection>(ps.getParameter<edm::InputTag>("uncleanBcCollection"));
0082   uncleanScCollection_ = consumes<SuperClusterCollection>(ps.getParameter<edm::InputTag>("uncleanScCollection"));
0083
0084   // the names of the products to be produced:
0085   //
0086   // the clean collection: this is as it was before, but labeled
0087   bcCollection_ = ps.getParameter<std::string>("bcCollection");
0088   scCollection_ = ps.getParameter<std::string>("scCollection");
0089   // the unclean only collection: SC unique to the unclean collection
0090   bcCollectionUncleanOnly_ = ps.getParameter<std::string>("bcCollectionUncleanOnly");
0091   scCollectionUncleanOnly_ = ps.getParameter<std::string>("scCollectionUncleanOnly");
0092   // the products:
0093   produces<reco::BasicClusterCollection>(bcCollection_);
0094   produces<reco::SuperClusterCollection>(scCollection_);
0095   produces<reco::BasicClusterCollection>(bcCollectionUncleanOnly_);
0096   produces<reco::SuperClusterCollection>(scCollectionUncleanOnly_);
0097 }
0098
0099 void UnifiedSCCollectionProducer::produce(edm::Event& evt, const edm::EventSetup& es) {
0100   edm::LogInfo("UnifiedSC") << ">>>>> Entering UnifiedSCCollectionProducer <<<<<";
0101   // get the input collections
0102   // __________________________________________________________________________
0103   //
0104   // cluster collections:
0105   edm::Handle<reco::BasicClusterCollection> pCleanBC;
0106   edm::Handle<reco::SuperClusterCollection> pCleanSC;
0107   //
0108   edm::Handle<reco::BasicClusterCollection> pUncleanBC;
0109   edm::Handle<reco::SuperClusterCollection> pUncleanSC;
0110
0111   evt.getByToken(cleanScCollection_, pCleanSC);
0112   evt.getByToken(cleanBcCollection_, pCleanBC);
0113   evt.getByToken(uncleanBcCollection_, pUncleanBC);
0114   evt.getByToken(uncleanScCollection_, pUncleanSC);
0115
0116   // the collections to be produced ___________________________________________
0117   reco::BasicClusterCollection basicClusters;
0118   reco::SuperClusterCollection superClusters;
0119   //
0120   reco::BasicClusterCollection basicClustersUncleanOnly;
0121   reco::SuperClusterCollection superClustersUncleanOnly;
0122   //
0123   // run over the uncleaned SC and check how many of them are matched to
0124   // the cleaned ones
0125   // if you find a matched one, then keep the info that it is matched
0126   //    along with which clean SC was matched + its basic clusters
0127   // if you find an unmatched one, keep the info and store its basic clusters
0128   //
0129   //
0130   int uncleanSize = pUncleanSC->size();
0131   int cleanSize = pCleanSC->size();
0132
0133   LogTrace("UnifiedSC") << "Size of Clean Collection: " << cleanSize << ", uncleanSize: " << uncleanSize;
0134
0135   // keep the indices
0136   std::vector<int> inUncleanOnlyInd;      // counting the unclean
0137   std::vector<int> inCleanInd;            // counting the unclean
0138   std::vector<int> inCleanOnlyInd;        // counting the clean
0139   std::vector<DetId> scUncleanSeedDetId;  // counting the unclean
0140   std::vector<DetId> scCleanSeedDetId;    // counting the clean
0141   // ontains the index of the SC that owns that BS
0142   // first basic cluster index, second: 0 for unclean and 1 for clean
0143   std::vector<std::pair<int, int> > basicClusterOwner;
0144   std::vector<std::pair<int, int> > basicClusterOwnerUncleanOnly;
0145   // if this basic cluster is a seed it is 1
0146   std::vector<int> uncleanBasicClusterIsSeed;
0147
0148   // loop over unclean SC _____________________________________________________
0149   for (int isc = 0; isc < uncleanSize; ++isc) {
0150     reco::SuperClusterRef unscRef(pUncleanSC, isc);
0151     const std::vector<std::pair<DetId, float> >& uhits = unscRef->hitsAndFractions();
0152     int uhitsSize = uhits.size();
0153     bool foundTheSame = false;
0154     for (int jsc = 0; jsc < cleanSize; ++jsc) {  // loop over the cleaned SC
0155       reco::SuperClusterRef cscRef(pCleanSC, jsc);
0156       const std::vector<std::pair<DetId, float> >& chits = cscRef->hitsAndFractions();
0157       int chitsSize = chits.size();
0158       foundTheSame = false;
0159       if (unscRef->seed()->seed() == cscRef->seed()->seed() && chitsSize == uhitsSize) {
0160         // if the clusters are exactly the same then because the clustering
0161         // algorithm works in a deterministic way, the order of the rechits
0162         // will be the same
0163         foundTheSame = true;
0164         for (int i = 0; i < chitsSize; ++i) {
0165           if (uhits[i].first != chits[i].first) {
0166             foundTheSame = false;
0167             break;
0168           }
0169         }
0170       }
0171       if (foundTheSame) {  // ok you have found it:
0172         // this supercluster belongs to both collections
0173         inUncleanOnlyInd.push_back(0);
0174         inCleanInd.push_back(jsc);  // keeps the index of the clean SC
0175         scUncleanSeedDetId.push_back(unscRef->seed()->seed());
0176         //
0177         // keep its basic clusters:
0178         for (reco::CaloCluster_iterator bciter = unscRef->clustersBegin(); bciter != unscRef->clustersEnd(); ++bciter) {
0179           // the basic clusters
0180           basicClusters.push_back(**bciter);
0181           // index of the unclean SC
0182           basicClusterOwner.push_back(std::make_pair(isc, 0));
0183         }
0184         break;  // break the loop over unclean sc
0185       }
0186     }
0187     if (not foundTheSame) {  // this SC is only in the unclean collection
0188       // mark it as unique in the uncleaned
0189       inUncleanOnlyInd.push_back(1);
0190       scUncleanSeedDetId.push_back(unscRef->seed()->seed());
0191       // keep all its basic clusters
0192       for (reco::CaloCluster_iterator bciter = unscRef->clustersBegin(); bciter != unscRef->clustersEnd(); ++bciter) {
0193         // the basic clusters
0194         basicClustersUncleanOnly.push_back(**bciter);
0195         basicClusterOwnerUncleanOnly.push_back(std::make_pair(isc, 0));
0196       }
0197     }
0198   }  // loop over the unclean SC _______________________________________________
0199   //
0200   int inCleanSize = inCleanInd.size();
0201   //
0202   // loop over the clean SC, check that are not in common with the unclean
0203   // ones and then store their SC as before ___________________________________
0204   for (int jsc = 0; jsc < cleanSize; ++jsc) {
0205     // check whether this index is already in the common collection
0206     bool takenAlready = false;
0207     for (int j = 0; j < inCleanSize; ++j) {
0208       if (jsc == inCleanInd[j]) {
0209         takenAlready = true;
0210         break;
0211       }
0212     }
0213     if (takenAlready) {
0214       inCleanOnlyInd.push_back(0);
0215       scCleanSeedDetId.push_back(DetId(0));
0216       continue;
0217     }
0218     inCleanOnlyInd.push_back(1);
0219     reco::SuperClusterRef cscRef(pCleanSC, jsc);
0220     scCleanSeedDetId.push_back(cscRef->seed()->seed());
0221     for (reco::CaloCluster_iterator bciter = cscRef->clustersBegin(); bciter != cscRef->clustersEnd(); ++bciter) {
0222       // the basic clusters
0223       basicClusters.push_back(**bciter);
0224       basicClusterOwner.push_back(std::make_pair(jsc, 1));
0225     }
0226   }  // end loop over clean SC _________________________________________________
0227   //
0228   //
0229
0230   // Final check: in the endcap BC may exist that are not associated to SC,
0231   // we need to recover them as well (e.g. multi5x5 algo)
0232   // This is should be optimized (SA, 20110621)
0233
0234   // loop on original clean BC collection and see if the BC is missing from the new one
0235   for (reco::BasicClusterCollection::const_iterator bc = pCleanBC->begin(); bc != pCleanBC->end(); ++bc) {
0236     bool foundTheSame = false;
0237     for (reco::BasicClusterCollection::const_iterator cleanonly_bc = basicClusters.begin();
0238          cleanonly_bc != basicClusters.end();
0239          ++cleanonly_bc) {
0240       const std::vector<std::pair<DetId, float> >& chits = bc->hitsAndFractions();
0241       int chitsSize = chits.size();
0242
0243       const std::vector<std::pair<DetId, float> >& uhits = cleanonly_bc->hitsAndFractions();
0244       int uhitsSize = uhits.size();
0245
0246       if (cleanonly_bc->seed() == bc->seed() && chitsSize == uhitsSize) {
0247         foundTheSame = true;
0248         for (int i = 0; i < chitsSize; ++i) {
0249           if (uhits[i].first != chits[i].first) {
0250             foundTheSame = false;
0251             break;
0252           }
0253         }
0254       }
0255
0256     }  // loop on new clean BC collection
0257
0258     // clean basic cluster is not associated to SC and does not belong to the
0259     // new collection, add it
0260     if (!foundTheSame) {
0261       basicClusters.push_back(*bc);
0262       LogTrace("UnifiedSC") << "found BC to add that was not associated to any SC";
0263     }
0264
0265   }  // loop on original clean BC collection
0266
0267   // at this point we have the basic cluster collection ready
0268   // Up to index   basicClusterOwner.size() we have the BC owned by a SC
0269   // The remaining are BCs not owned by a SC
0270
0271   int bcSize = (int)basicClusterOwner.size();
0272   int bcSizeUncleanOnly = (int)basicClustersUncleanOnly.size();
0273
0274   LogTrace("UnifiedSC") << "Found cleaned SC: " << cleanSize << " uncleaned SC: " << uncleanSize;
0275   //
0276   // export the clusters to the event from the clean clusters
0277   auto basicClusters_p = std::make_unique<reco::BasicClusterCollection>();
0278   basicClusters_p->assign(basicClusters.begin(), basicClusters.end());
0279   edm::OrphanHandle<reco::BasicClusterCollection> bccHandle = evt.put(std::move(basicClusters_p), bcCollection_);
0280   if (!(bccHandle.isValid())) {
0281     edm::LogWarning("MissingInput") << "could not handle the new BasicClusters!";
0282     return;
0283   }
0284
0285   LogTrace("UnifiedSC") << "Got the BasicClusters from the event again";
0286   //
0287   // export the clusters to the event: from the unclean only clusters
0288   auto basicClustersUncleanOnly_p = std::make_unique<reco::BasicClusterCollection>();
0289   basicClustersUncleanOnly_p->assign(basicClustersUncleanOnly.begin(), basicClustersUncleanOnly.end());
0290   edm::OrphanHandle<reco::BasicClusterCollection> bccHandleUncleanOnly =
0291       evt.put(std::move(basicClustersUncleanOnly_p), bcCollectionUncleanOnly_);
0292   if (!(bccHandleUncleanOnly.isValid())) {
0293     edm::LogWarning("MissingInput") << "could not handle the new BasicClusters (Unclean Only)!";
0294     return;
0295   }
0296   LogTrace("UnifiedSC") << "Got the BasicClusters from the event again  (Unclean Only)";
0297   //
0298
0299   // now we can build the SC collection
0300   //
0301   // start again from the unclean collection
0302   // all the unclean SC will become members of the new collection
0303   // with different algoIDs ___________________________________________________
0304   for (int isc = 0; isc < uncleanSize; ++isc) {
0305     reco::CaloClusterPtrVector clusterPtrVector;
0306     // the seed is the basic cluster with the highest energy
0307     reco::CaloClusterPtr seed;
0308     if (inUncleanOnlyInd[isc] == 1) {  // unclean SC Unique in Unclean
0309       for (int jbc = 0; jbc < bcSizeUncleanOnly; ++jbc) {
0310         std::pair<int, int> theBcOwner = basicClusterOwnerUncleanOnly[jbc];
0311         if (theBcOwner.first == isc && theBcOwner.second == 0) {
0312           reco::CaloClusterPtr currentClu = reco::CaloClusterPtr(bccHandleUncleanOnly, jbc);
0313           clusterPtrVector.push_back(currentClu);
0314           if (scUncleanSeedDetId[isc] == currentClu->seed()) {
0315             seed = currentClu;
0316           }
0317         }
0318       }
0319
0320     } else {  // unclean SC common in clean and unclean
0321       for (int jbc = 0; jbc < bcSize; ++jbc) {
0322         std::pair<int, int> theBcOwner = basicClusterOwner[jbc];
0323         if (theBcOwner.first == isc && theBcOwner.second == 0) {
0324           reco::CaloClusterPtr currentClu = reco::CaloClusterPtr(bccHandle, jbc);
0325           clusterPtrVector.push_back(currentClu);
0326           if (scUncleanSeedDetId[isc] == currentClu->seed()) {
0327             seed = currentClu;
0328           }
0329         }
0330       }
0331     }
0332     //std::cout << "before getting the uncl" << std::endl;
0333     reco::SuperClusterRef unscRef(pUncleanSC, isc);
0334     reco::SuperCluster newSC(unscRef->energy(), unscRef->position(), seed, clusterPtrVector);
0335     // now set the algoID for this SC again
0336     if (inUncleanOnlyInd[isc] == 1) {
0337       // set up the quality to unclean only .............
0338       newSC.setFlags(reco::CaloCluster::uncleanOnly);
0339       superClustersUncleanOnly.push_back(newSC);
0340     } else {
0341       // set up the quality to common  .............
0342       newSC.setFlags(reco::CaloCluster::common);
0343       superClusters.push_back(newSC);
0344     }
0345     // now you can store your SC
0346
0347   }  // end loop over unclean SC _______________________________________________
0348   //  flags numbering scheme
0349   //  flags =   0 = cleanedOnly     cluster is only in the cleaned collection
0350   //  flags = 100 = common          cluster is common in both collections
0351   //  flags = 200 = uncleanedOnly   cluster is only in the uncleaned collection
0352
0353   // now loop over the clean SC and do the same but now you have to avoid the
0354   // the duplicated ones ______________________________________________________
0355   for (int jsc = 0; jsc < cleanSize; ++jsc) {
0356     //std::cout << "working in cl #" << jsc << std::endl;
0357     // check that the SC is not in the unclean collection
0358     if (inCleanOnlyInd[jsc] == 0)
0359       continue;
0360     reco::CaloClusterPtrVector clusterPtrVector;
0361     // the seed is the basic cluster with the highest energy
0362     reco::CaloClusterPtr seed;
0363     for (int jbc = 0; jbc < bcSize; ++jbc) {
0364       std::pair<int, int> theBcOwner = basicClusterOwner[jbc];
0365       if (theBcOwner.first == jsc && theBcOwner.second == 1) {
0366         reco::CaloClusterPtr currentClu = reco::CaloClusterPtr(bccHandle, jbc);
0367         clusterPtrVector.push_back(currentClu);
0368         if (scCleanSeedDetId[jsc] == currentClu->seed()) {
0369           seed = currentClu;
0370         }
0371       }
0372     }
0373     reco::SuperClusterRef cscRef(pCleanSC, jsc);
0374     reco::SuperCluster newSC(cscRef->energy(), cscRef->position(), seed, clusterPtrVector);
0375     newSC.setFlags(reco::CaloCluster::cleanOnly);
0376
0377     // add it to the collection:
0378     superClusters.push_back(newSC);
0379
0380   }  // end loop over clean SC _________________________________________________
0381
0382   LogTrace("UnifiedSC") << "New SC collection was created";
0383
0384   auto superClusters_p = std::make_unique<reco::SuperClusterCollection>();
0385   superClusters_p->assign(superClusters.begin(), superClusters.end());
0386
0387   evt.put(std::move(superClusters_p), scCollection_);
0388
0389   LogTrace("UnifiedSC") << "Clusters (Basic/Super) added to the Event! :-)";
0390
0391   auto superClustersUncleanOnly_p = std::make_unique<reco::SuperClusterCollection>();
0392   superClustersUncleanOnly_p->assign(superClustersUncleanOnly.begin(), superClustersUncleanOnly.end());
0393
0394   evt.put(std::move(superClustersUncleanOnly_p), scCollectionUncleanOnly_);
0395
0396   // ----- debugging ----------
0397   // print the new collection SC quantities
0398
0399   // print out the clean collection SC
0400   LogTrace("UnifiedSC") << "Clean Collection SC ";
0401   for (int i = 0; i < cleanSize; ++i) {
0402     reco::SuperClusterRef cscRef(pCleanSC, i);
0403     LogTrace("UnifiedSC") << " >>> clean    #" << i << "; Energy: " << cscRef->energy() << " eta: " << cscRef->eta()
0404                           << " sc seed detid: " << cscRef->seed()->seed().rawId() << std::endl;
0405   }
0406   // the unclean SC
0407   LogTrace("UnifiedSC") << "Unclean Collection SC ";
0408   for (int i = 0; i < uncleanSize; ++i) {
0409     reco::SuperClusterRef uscRef(pUncleanSC, i);
0410     LogTrace("UnifiedSC") << " >>> unclean  #" << i << "; Energy: " << uscRef->energy() << " eta: " << uscRef->eta()
0411                           << " sc seed detid: " << uscRef->seed()->seed().rawId();
0412   }
0413   // the new collection
0414   LogTrace("UnifiedSC") << "The new SC clean collection with size " << superClusters.size() << std::endl;
0415
0416   int new_unclean = 0, new_clean = 0;
0417   for (int i = 0; i < (int)superClusters.size(); ++i) {
0418     const reco::SuperCluster& nsc = superClusters[i];
0419     LogTrace("UnifiedSC") << "SC was got" << std::endl
0420                           << " ---> energy: " << nsc.energy() << std::endl
0421                           << " ---> eta: " << nsc.eta() << std::endl
0422                           << " ---> inClean: " << nsc.isInClean() << std::endl
0423                           << " ---> id: " << nsc.seed()->seed().rawId() << std::endl
0424                           << " >>> newSC    #" << i << "; Energy: " << nsc.energy() << " eta: " << nsc.eta()
0425                           << " isClean=" << nsc.isInClean() << " isUnclean=" << nsc.isInUnclean()
0426                           << " sc seed detid: " << nsc.seed()->seed().rawId();
0427
0428     if (nsc.isInUnclean())
0429       ++new_unclean;
0430     if (nsc.isInClean())
0431       ++new_clean;
0432   }
0433   LogTrace("UnifiedSC") << "The new SC unclean only collection with size " << superClustersUncleanOnly.size();
0434   for (int i = 0; i < (int)superClustersUncleanOnly.size(); ++i) {
0435     const reco::SuperCluster nsc = superClustersUncleanOnly[i];
0436     LogTrace("UnifiedSC") << " >>> newSC    #" << i << "; Energy: " << nsc.energy() << " eta: " << nsc.eta()
0437                           << " isClean=" << nsc.isInClean() << " isUnclean=" << nsc.isInUnclean()
0438                           << " sc seed detid: " << nsc.seed()->seed().rawId();
0439     if (nsc.isInUnclean())
0440       ++new_unclean;
0441     if (nsc.isInClean())
0442       ++new_clean;
0443   }
0444   if ((new_unclean != uncleanSize) || (new_clean != cleanSize)) {
0445     LogTrace("UnifiedSC") << ">>>>!!!!!! MISMATCH: new unclean/ old unclean= " << new_unclean << " / " << uncleanSize
0446                           << ", new clean/ old clean" << new_clean << " / " << cleanSize;
0447   }
0448 }