JetMCAlgos/plugins/CandOneToOneDeltaRMatcher.cc

0001 /* \class CandOneToOneDeltaRMatcher
0002  *
0003  * Producer for simple match map
0004  * to match two collections of candidate
0005  * with one-to-One matching
0006  * minimizing Sum(DeltaR)
0007  *
0008  */
0009
0010 #include "FWCore/Framework/interface/global/EDProducer.h"
0011 #include "FWCore/ParameterSet/interface/ParameterSetfwd.h"
0012
0013 #include "DataFormats/Candidate/interface/Candidate.h"
0014
0015 #include <vector>
0016 #include <iostream>
0017
0018 class CandOneToOneDeltaRMatcher : public edm::global::EDProducer<> {
0019 public:
0020   explicit CandOneToOneDeltaRMatcher(const edm::ParameterSet&);
0021
0022   using AllDist = std::vector<std::vector<float>>;
0023
0024   enum class Algo { kBruteForce, kSwitchMode, kMixMode };
0025
0026 private:
0027   void produce(edm::StreamID, edm::Event&, const edm::EventSetup&) const override;
0028   double length(const std::vector<int>&, const AllDist&) const;
0029   std::vector<int> AlgoBruteForce(int, int, const AllDist&) const;
0030   std::vector<int> AlgoSwitchMethod(int, int, AllDist&) const;
0031   static Algo algo(std::string const&);
0032   static const char* algoName(Algo);
0033
0034   const edm::EDGetTokenT<reco::CandidateView> sourceToken_;
0035   const edm::EDGetTokenT<reco::CandidateView> matchedToken_;
0036   const Algo algoMethod_;
0037 };
0038
0039 #include "PhysicsTools/JetMCUtils/interface/combination.h"
0040
0041 #include "FWCore/Framework/interface/ESHandle.h"
0042 #include "FWCore/Framework/interface/Event.h"
0043 #include "FWCore/Framework/interface/EventSetup.h"
0044 #include "FWCore/Utilities/interface/InputTag.h"
0045 #include "FWCore/Utilities/interface/EDMException.h"
0046 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0047 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0048
0049 #include "DataFormats/Common/interface/Handle.h"
0050 #include "DataFormats/Candidate/interface/LeafCandidate.h"
0051 #include "DataFormats/Candidate/interface/CandMatchMap.h"
0052 #include "DataFormats/Candidate/interface/CandidateFwd.h"
0053
0054 #include <Math/VectorUtil.h>
0055 #include <TMath.h>
0056
0057 using namespace edm;
0058 using namespace std;
0059 using namespace reco;
0060 using namespace ROOT::Math::VectorUtil;
0061 using namespace stdcomb;
0062
0063 CandOneToOneDeltaRMatcher::Algo CandOneToOneDeltaRMatcher::algo(const std::string& algoName) {
0064   if (algoName == "BruteForce") {
0065     return Algo::kBruteForce;
0066   } else if (algoName == "SwitchMode") {
0067     return Algo::kSwitchMode;
0068   } else if (algoName == "MixMode") {
0069     return Algo::kMixMode;
0070   } else {
0071     throw cms::Exception("OneToOne Constructor") << "wrong matching method in ParameterSet";
0072   }
0073 }
0074
0075 const char* CandOneToOneDeltaRMatcher::algoName(Algo iAlgo) {
0076   switch (iAlgo) {
0077     case Algo::kBruteForce:
0078       return "BruteForce";
0079     case Algo::kSwitchMode:
0080       return "SwitchMode";
0081     case Algo::kMixMode:
0082       return "MixMode";
0083   }
0084   //can not get here
0085   return "";
0086 }
0087
0088 CandOneToOneDeltaRMatcher::CandOneToOneDeltaRMatcher(const ParameterSet& cfg)
0089     : sourceToken_(consumes<CandidateView>(cfg.getParameter<InputTag>("src"))),
0090       matchedToken_(consumes<CandidateView>(cfg.getParameter<InputTag>("matched"))),
0091       algoMethod_(algo(cfg.getParameter<string>("algoMethod"))) {
0092   produces<CandViewMatchMap>("src2mtc");
0093   produces<CandViewMatchMap>("mtc2src");
0094 }
0095
0096 void CandOneToOneDeltaRMatcher::produce(StreamID, Event& evt, const EventSetup& es) const {
0097   Handle<CandidateView> source;
0098   Handle<CandidateView> matched;
0099   evt.getByToken(sourceToken_, source);
0100   evt.getByToken(matchedToken_, matched);
0101
0102   edm::LogVerbatim("CandOneToOneDeltaRMatcher") << "======== Source Collection =======";
0103   for (CandidateView::const_iterator c = source->begin(); c != source->end(); ++c) {
0104     edm::LogVerbatim("CandOneToOneDeltaRMatcher")
0105         << " pt source  " << c->pt() << " " << c->eta() << " " << c->phi() << endl;
0106   }
0107   edm::LogVerbatim("CandOneToOneDeltaRMatcher") << "======== Matched Collection =======";
0108   for (CandidateView::const_iterator c = matched->begin(); c != matched->end(); ++c) {
0109     edm::LogVerbatim("CandOneToOneDeltaRMatcher")
0110         << " pt source  " << c->pt() << " " << c->eta() << " " << c->phi() << endl;
0111   }
0112
0113   const int nSrc = source->size();
0114   const int nMtc = matched->size();
0115
0116   const int nMin = min(source->size(), matched->size());
0117   const int nMax = max(source->size(), matched->size());
0118   if (nMin < 1)
0119     return;
0120
0121   std::vector<std::vector<float>> allDist;
0122
0123   if (nSrc <= nMtc) {
0124     allDist.reserve(source->size());
0125     for (CandidateView::const_iterator iSr = source->begin(); iSr != source->end(); iSr++) {
0126       vector<float> tempAllDist;
0127       tempAllDist.reserve(matched->size());
0128       for (CandidateView::const_iterator iMt = matched->begin(); iMt != matched->end(); iMt++) {
0129         tempAllDist.push_back(DeltaR(iSr->p4(), iMt->p4()));
0130       }
0131       allDist.emplace_back(std::move(tempAllDist));
0132     }
0133   } else {
0134     allDist.reserve(matched->size());
0135     for (CandidateView::const_iterator iMt = matched->begin(); iMt != matched->end(); iMt++) {
0136       vector<float> tempAllDist;
0137       tempAllDist.reserve(source->size());
0138       for (CandidateView::const_iterator iSr = source->begin(); iSr != source->end(); iSr++) {
0139         tempAllDist.push_back(DeltaR(iSr->p4(), iMt->p4()));
0140       }
0141       allDist.emplace_back(std::move(tempAllDist));
0142     }
0143   }
0144
0145   /*
0146   edm::LogVerbatim("CandOneToOneDeltaRMatcher") << "======== The DeltaR Matrix =======";
0147   for(int m0=0; m0<nMin; m0++) {
0148     //    for(int m1=0; m1<nMax; m1++) {
0149       edm::LogVerbatim("CandOneToOneDeltaRMatcher") << setprecision(2) << fixed << (m1 AllDist[m0][m1] ;
0150     //}
0151     edm::LogVerbatim("CandOneToOneDeltaRMatcher") << "\n";
0152   }
0153   */
0154
0155   // Loop size if Brute Force
0156   int nLoopToDo = (int)(TMath::Factorial(nMax) / TMath::Factorial(nMax - nMin));
0157   edm::LogVerbatim("CandOneToOneDeltaRMatcher") << "nLoop:" << nLoopToDo << endl;
0158   edm::LogVerbatim("CandOneToOneDeltaRMatcher") << "Choosen Algo is:" << algoName(algoMethod_);
0159   vector<int> bestCB;
0160
0161   // Algo is Brute Force
0162   if (algoMethod_ == Algo::kBruteForce) {
0163     bestCB = AlgoBruteForce(nMin, nMax, allDist);
0164
0165     // Algo is Switch Method
0166   } else if (algoMethod_ == Algo::kSwitchMode) {
0167     bestCB = AlgoSwitchMethod(nMin, nMax, allDist);
0168
0169     // Algo is Brute Force if nLoop < 10000
0170   } else if (algoMethod_ == Algo::kMixMode) {
0171     if (nLoopToDo < 10000) {
0172       bestCB = AlgoBruteForce(nMin, nMax, allDist);
0173     } else {
0174       bestCB = AlgoSwitchMethod(nMin, nMax, allDist);
0175     }
0176   }
0177
0178   for (int i1 = 0; i1 < nMin; i1++)
0179     edm::LogVerbatim("CandOneToOneDeltaRMatcher")
0180         << "min: " << i1 << " " << bestCB[i1] << " " << allDist[i1][bestCB[i1]];
0181
0182   /*
0183   auto matchMapSrMt = std::make_unique<CandViewMatchMap>(CandViewMatchMap::ref_type( CandidateRefProd( source  ),
0184                                                                                              CandidateRefProd( matched ) ) );
0185   auto matchMapMtSr = std::make_unique<CandViewMatchMap>(CandViewMatchMap::ref_type( CandidateRefProd( matched ),
0186                                                                                              CandidateRefProd( source ) ) );
0187 */
0188
0189   auto matchMapSrMt = std::make_unique<CandViewMatchMap>();
0190   auto matchMapMtSr = std::make_unique<CandViewMatchMap>();
0191
0192   for (int c = 0; c != nMin; c++) {
0193     if (source->size() <= matched->size()) {
0194       matchMapSrMt->insert(source->refAt(c), matched->refAt(bestCB[c]));
0195       matchMapMtSr->insert(matched->refAt(bestCB[c]), source->refAt(c));
0196     } else {
0197       matchMapSrMt->insert(source->refAt(bestCB[c]), matched->refAt(c));
0198       matchMapMtSr->insert(matched->refAt(c), source->refAt(bestCB[c]));
0199     }
0200   }
0201
0202   /*
0203   for( int c = 0; c != nMin; c ++ ) {
0204     if( source->size() <= matched->size() ) {
0205       matchMapSrMt->insert( CandidateRef( source,  c         ), CandidateRef( matched, bestCB[c] ) );
0206       matchMapMtSr->insert( CandidateRef( matched, bestCB[c] ), CandidateRef( source, c          ) );
0207     } else {
0208       matchMapSrMt->insert( CandidateRef( source,  bestCB[c] ), CandidateRef( matched, c         ) );
0209       matchMapMtSr->insert( CandidateRef( matched, c         ), CandidateRef( source,  bestCB[c] ) );
0210     }
0211   }
0212 */
0213   evt.put(std::move(matchMapSrMt), "src2mtc");
0214   evt.put(std::move(matchMapMtSr), "mtc2src");
0215 }
0216
0217 double CandOneToOneDeltaRMatcher::length(const vector<int>& best, const AllDist& allDist) const {
0218   double myLength = 0;
0219   int row = 0;
0220   for (vector<int>::const_iterator it = best.begin(); it != best.end(); it++) {
0221     myLength += allDist[row][*it];
0222     row++;
0223   }
0224   return myLength;
0225 }
0226
0227 // this is the Brute Force Algorithm
0228 // All the possible combination are checked
0229 // The best one is always found
0230 // Be carefull when you have high values for nMin and nMax --> the combinatorial could explode!
0231 // Sum(DeltaR) is minimized -->
0232 // 0.1 - 0.2 - 1.0 - 1.5 is lower than
0233 // 0.1 - 0.2 - 0.3 - 3.0
0234 // Which one do you prefer? --> BruteForce select always the first
0235
0236 vector<int> CandOneToOneDeltaRMatcher::AlgoBruteForce(int nMin, int nMax, const AllDist& allDist) const {
0237   vector<int> ca;
0238   vector<int> cb;
0239   vector<int> bestCB;
0240   float totalDeltaR = 0;
0241   float BestTotalDeltaR = 1000;
0242
0243   ca.reserve(nMax);
0244   for (int i1 = 0; i1 < nMax; i1++)
0245     ca.push_back(i1);
0246   cb.reserve(nMin);
0247   for (int i1 = 0; i1 < nMin; i1++)
0248     cb.push_back(i1);
0249
0250   do {
0251     //do your processing on the new combination here
0252     for (int cnt = 0; cnt < TMath::Factorial(nMin); cnt++) {
0253       totalDeltaR = length(cb, allDist);
0254       if (totalDeltaR < BestTotalDeltaR) {
0255         BestTotalDeltaR = totalDeltaR;
0256         bestCB = cb;
0257       }
0258       next_permutation(cb.begin(), cb.end());
0259     }
0260   } while (next_combination(ca.begin(), ca.end(), cb.begin(), cb.end()));
0261
0262   return bestCB;
0263 }
0264
0265 // This method (Developed originally by Daniele Benedetti) check for the best combination
0266 // choosing the minimum DeltaR for each line in AllDist matrix
0267 // If no repeated row is found: ie (line,col)=(1,3) and (2,3) --> same as BruteForce
0268 // If repetition --> set the higher DeltaR between  the 2 repetition to 1000 and re-check best combination
0269 // Iterate until no repetition
0270 // No guaranted minimum for Sum(DeltaR)
0271 // If you have:
0272 // 0.1 - 0.2 - 1.0 - 1.5 is lower than
0273 // 0.1 - 0.2 - 0.3 - 3.0
0274 // SwitchMethod normally select the second solution
0275
0276 vector<int> CandOneToOneDeltaRMatcher::AlgoSwitchMethod(int nMin, int nMax, AllDist& allDist) const {
0277   vector<int> bestCB;
0278   for (int i1 = 0; i1 < nMin; i1++) {
0279     int minInd = 0;
0280     for (int i2 = 1; i2 < nMax; i2++)
0281       if (allDist[i1][i2] < allDist[i1][minInd])
0282         minInd = i2;
0283     bestCB.push_back(minInd);
0284   }
0285
0286   bool inside = true;
0287   while (inside) {
0288     inside = false;
0289     for (int i1 = 0; i1 < nMin; i1++) {
0290       for (int i2 = i1 + 1; i2 < nMin; i2++) {
0291         if (bestCB[i1] == bestCB[i2]) {
0292           inside = true;
0293           if (allDist[i1][(bestCB[i1])] <= allDist[i2][(bestCB[i2])]) {
0294             allDist[i2][(bestCB[i2])] = 1000;
0295             int minInd = 0;
0296             for (int i3 = 1; i3 < nMax; i3++)
0297               if (allDist[i2][i3] < allDist[i2][minInd])
0298                 minInd = i3;
0299             bestCB[i2] = minInd;
0300           } else {
0301             allDist[i1][(bestCB[i1])] = 1000;
0302             int minInd = 0;
0303             for (int i3 = 1; i3 < nMax; i3++)
0304               if (allDist[i1][i3] < allDist[i1][minInd])
0305                 minInd = i3;
0306             bestCB[i1] = minInd;
0307           }
0308         }  // End if
0309       }
0310     }
0311   }  // End while
0312
0313   return bestCB;
0314 }
0315
0316 #include "FWCore/PluginManager/interface/ModuleDef.h"
0317 #include "FWCore/Framework/interface/MakerMacros.h"
0318
0319 DEFINE_FWK_MODULE(CandOneToOneDeltaRMatcher);