Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-06-08 05:47:15

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   reco::BasicClusterCollection basicClustersProd = *bccHandle;
0285 
0286   LogTrace("UnifiedSC") << "Got the BasicClusters from the event again";
0287   //
0288   // export the clusters to the event: from the unclean only clusters
0289   auto basicClustersUncleanOnly_p = std::make_unique<reco::BasicClusterCollection>();
0290   basicClustersUncleanOnly_p->assign(basicClustersUncleanOnly.begin(), basicClustersUncleanOnly.end());
0291   edm::OrphanHandle<reco::BasicClusterCollection> bccHandleUncleanOnly =
0292       evt.put(std::move(basicClustersUncleanOnly_p), bcCollectionUncleanOnly_);
0293   if (!(bccHandleUncleanOnly.isValid())) {
0294     edm::LogWarning("MissingInput") << "could not handle the new BasicClusters (Unclean Only)!";
0295     return;
0296   }
0297   reco::BasicClusterCollection basicClustersUncleanOnlyProd = *bccHandleUncleanOnly;
0298   LogTrace("UnifiedSC") << "Got the BasicClusters from the event again  (Unclean Only)";
0299   //
0300 
0301   // now we can build the SC collection
0302   //
0303   // start again from the unclean collection
0304   // all the unclean SC will become members of the new collection
0305   // with different algoIDs ___________________________________________________
0306   for (int isc = 0; isc < uncleanSize; ++isc) {
0307     reco::CaloClusterPtrVector clusterPtrVector;
0308     // the seed is the basic cluster with the highest energy
0309     reco::CaloClusterPtr seed;
0310     if (inUncleanOnlyInd[isc] == 1) {  // unclean SC Unique in Unclean
0311       for (int jbc = 0; jbc < bcSizeUncleanOnly; ++jbc) {
0312         std::pair<int, int> theBcOwner = basicClusterOwnerUncleanOnly[jbc];
0313         if (theBcOwner.first == isc && theBcOwner.second == 0) {
0314           reco::CaloClusterPtr currentClu = reco::CaloClusterPtr(bccHandleUncleanOnly, jbc);
0315           clusterPtrVector.push_back(currentClu);
0316           if (scUncleanSeedDetId[isc] == currentClu->seed()) {
0317             seed = currentClu;
0318           }
0319         }
0320       }
0321 
0322     } else {  // unclean SC common in clean and unclean
0323       for (int jbc = 0; jbc < bcSize; ++jbc) {
0324         std::pair<int, int> theBcOwner = basicClusterOwner[jbc];
0325         if (theBcOwner.first == isc && theBcOwner.second == 0) {
0326           reco::CaloClusterPtr currentClu = reco::CaloClusterPtr(bccHandle, jbc);
0327           clusterPtrVector.push_back(currentClu);
0328           if (scUncleanSeedDetId[isc] == currentClu->seed()) {
0329             seed = currentClu;
0330           }
0331         }
0332       }
0333     }
0334     //std::cout << "before getting the uncl" << std::endl;
0335     reco::SuperClusterRef unscRef(pUncleanSC, isc);
0336     reco::SuperCluster newSC(unscRef->energy(), unscRef->position(), seed, clusterPtrVector);
0337     // now set the algoID for this SC again
0338     if (inUncleanOnlyInd[isc] == 1) {
0339       // set up the quality to unclean only .............
0340       newSC.setFlags(reco::CaloCluster::uncleanOnly);
0341       superClustersUncleanOnly.push_back(newSC);
0342     } else {
0343       // set up the quality to common  .............
0344       newSC.setFlags(reco::CaloCluster::common);
0345       superClusters.push_back(newSC);
0346     }
0347     // now you can store your SC
0348 
0349   }  // end loop over unclean SC _______________________________________________
0350   //  flags numbering scheme
0351   //  flags =   0 = cleanedOnly     cluster is only in the cleaned collection
0352   //  flags = 100 = common          cluster is common in both collections
0353   //  flags = 200 = uncleanedOnly   cluster is only in the uncleaned collection
0354 
0355   // now loop over the clean SC and do the same but now you have to avoid the
0356   // the duplicated ones ______________________________________________________
0357   for (int jsc = 0; jsc < cleanSize; ++jsc) {
0358     //std::cout << "working in cl #" << jsc << std::endl;
0359     // check that the SC is not in the unclean collection
0360     if (inCleanOnlyInd[jsc] == 0)
0361       continue;
0362     reco::CaloClusterPtrVector clusterPtrVector;
0363     // the seed is the basic cluster with the highest energy
0364     reco::CaloClusterPtr seed;
0365     for (int jbc = 0; jbc < bcSize; ++jbc) {
0366       std::pair<int, int> theBcOwner = basicClusterOwner[jbc];
0367       if (theBcOwner.first == jsc && theBcOwner.second == 1) {
0368         reco::CaloClusterPtr currentClu = reco::CaloClusterPtr(bccHandle, jbc);
0369         clusterPtrVector.push_back(currentClu);
0370         if (scCleanSeedDetId[jsc] == currentClu->seed()) {
0371           seed = currentClu;
0372         }
0373       }
0374     }
0375     reco::SuperClusterRef cscRef(pCleanSC, jsc);
0376     reco::SuperCluster newSC(cscRef->energy(), cscRef->position(), seed, clusterPtrVector);
0377     newSC.setFlags(reco::CaloCluster::cleanOnly);
0378 
0379     // add it to the collection:
0380     superClusters.push_back(newSC);
0381 
0382   }  // end loop over clean SC _________________________________________________
0383 
0384   LogTrace("UnifiedSC") << "New SC collection was created";
0385 
0386   auto superClusters_p = std::make_unique<reco::SuperClusterCollection>();
0387   superClusters_p->assign(superClusters.begin(), superClusters.end());
0388 
0389   evt.put(std::move(superClusters_p), scCollection_);
0390 
0391   LogTrace("UnifiedSC") << "Clusters (Basic/Super) added to the Event! :-)";
0392 
0393   auto superClustersUncleanOnly_p = std::make_unique<reco::SuperClusterCollection>();
0394   superClustersUncleanOnly_p->assign(superClustersUncleanOnly.begin(), superClustersUncleanOnly.end());
0395 
0396   evt.put(std::move(superClustersUncleanOnly_p), scCollectionUncleanOnly_);
0397 
0398   // ----- debugging ----------
0399   // print the new collection SC quantities
0400 
0401   // print out the clean collection SC
0402   LogTrace("UnifiedSC") << "Clean Collection SC ";
0403   for (int i = 0; i < cleanSize; ++i) {
0404     reco::SuperClusterRef cscRef(pCleanSC, i);
0405     LogTrace("UnifiedSC") << " >>> clean    #" << i << "; Energy: " << cscRef->energy() << " eta: " << cscRef->eta()
0406                           << " sc seed detid: " << cscRef->seed()->seed().rawId() << std::endl;
0407   }
0408   // the unclean SC
0409   LogTrace("UnifiedSC") << "Unclean Collection SC ";
0410   for (int i = 0; i < uncleanSize; ++i) {
0411     reco::SuperClusterRef uscRef(pUncleanSC, i);
0412     LogTrace("UnifiedSC") << " >>> unclean  #" << i << "; Energy: " << uscRef->energy() << " eta: " << uscRef->eta()
0413                           << " sc seed detid: " << uscRef->seed()->seed().rawId();
0414   }
0415   // the new collection
0416   LogTrace("UnifiedSC") << "The new SC clean collection with size " << superClusters.size() << std::endl;
0417 
0418   int new_unclean = 0, new_clean = 0;
0419   for (int i = 0; i < (int)superClusters.size(); ++i) {
0420     const reco::SuperCluster& nsc = superClusters[i];
0421     LogTrace("UnifiedSC") << "SC was got" << std::endl
0422                           << " ---> energy: " << nsc.energy() << std::endl
0423                           << " ---> eta: " << nsc.eta() << std::endl
0424                           << " ---> inClean: " << nsc.isInClean() << std::endl
0425                           << " ---> id: " << nsc.seed()->seed().rawId() << std::endl
0426                           << " >>> newSC    #" << i << "; Energy: " << nsc.energy() << " eta: " << nsc.eta()
0427                           << " isClean=" << nsc.isInClean() << " isUnclean=" << nsc.isInUnclean()
0428                           << " sc seed detid: " << nsc.seed()->seed().rawId();
0429 
0430     if (nsc.isInUnclean())
0431       ++new_unclean;
0432     if (nsc.isInClean())
0433       ++new_clean;
0434   }
0435   LogTrace("UnifiedSC") << "The new SC unclean only collection with size " << superClustersUncleanOnly.size();
0436   for (int i = 0; i < (int)superClustersUncleanOnly.size(); ++i) {
0437     const reco::SuperCluster nsc = superClustersUncleanOnly[i];
0438     LogTrace("UnifiedSC") << " >>> newSC    #" << i << "; Energy: " << nsc.energy() << " eta: " << nsc.eta()
0439                           << " isClean=" << nsc.isInClean() << " isUnclean=" << nsc.isInUnclean()
0440                           << " sc seed detid: " << nsc.seed()->seed().rawId();
0441     if (nsc.isInUnclean())
0442       ++new_unclean;
0443     if (nsc.isInClean())
0444       ++new_clean;
0445   }
0446   if ((new_unclean != uncleanSize) || (new_clean != cleanSize)) {
0447     LogTrace("UnifiedSC") << ">>>>!!!!!! MISMATCH: new unclean/ old unclean= " << new_unclean << " / " << uncleanSize
0448                           << ", new clean/ old clean" << new_clean << " / " << cleanSize;
0449   }
0450 }