Reference Manual

Go to the documentation of this file.
00001 /* \class CandOneToOneDeltaRMatcher
00002  *
00003  * Producer for simple match map
00004  * to match two collections of candidate
00005  * with one-to-One matching 
00006  * minimizing Sum(DeltaR)
00007  *
00008  */
00009 
00010 #include "FWCore/Framework/interface/EDProducer.h"
00011 #include "FWCore/ParameterSet/interface/ParameterSetfwd.h"
00012 #include "FWCore/Utilities/interface/InputTag.h"
00013 
00014 #include<vector>
00015 #include<iostream>
00016 
00017 class CandOneToOneDeltaRMatcher : public edm::EDProducer {
00018  public:
00019   CandOneToOneDeltaRMatcher( const edm::ParameterSet & );
00020   ~CandOneToOneDeltaRMatcher();
00021  private:
00022   void produce( edm::Event&, const edm::EventSetup& );
00023   double lenght( const std::vector<int>& );
00024   std::vector<int> AlgoBruteForce(int, int);
00025   std::vector<int> AlgoSwitchMethod(int, int);
00026   
00027   edm::InputTag source_;
00028   edm::InputTag matched_;
00029   std::vector < std::vector<float> > AllDist;
00030   std::string algoMethod_;
00031 
00032 };
00033 
00034 #include "PhysicsTools/JetMCUtils/interface/combination.h"
00035 
00036 #include "FWCore/Framework/interface/ESHandle.h"
00037 #include "FWCore/Framework/interface/Event.h"
00038 #include "FWCore/Framework/interface/EventSetup.h"
00039 #include "FWCore/Utilities/interface/EDMException.h"
00040 #include "FWCore/ParameterSet/interface/ParameterSet.h"
00041 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00042 
00043 #include "DataFormats/Common/interface/Handle.h"
00044 #include "DataFormats/Candidate/interface/Candidate.h"
00045 #include "DataFormats/Candidate/interface/LeafCandidate.h"
00046 #include "DataFormats/Candidate/interface/CandMatchMap.h"
00047 #include "DataFormats/Candidate/interface/CandidateFwd.h"
00048 
00049 
00050 #include <Math/VectorUtil.h>
00051 #include <TMath.h>
00052 
00053 using namespace edm;
00054 using namespace std;
00055 using namespace reco;
00056 using namespace ROOT::Math::VectorUtil;
00057 using namespace stdcomb;
00058 
00059 CandOneToOneDeltaRMatcher::CandOneToOneDeltaRMatcher( const ParameterSet & cfg ) :
00060   source_( cfg.getParameter<InputTag>( "src" ) ),
00061   matched_( cfg.getParameter<InputTag>( "matched" ) ),
00062   algoMethod_( cfg.getParameter<string>( "algoMethod" ) ) {
00063   produces<CandViewMatchMap>("src2mtc");
00064   produces<CandViewMatchMap>("mtc2src");
00065 }
00066 
00067 CandOneToOneDeltaRMatcher::~CandOneToOneDeltaRMatcher() {
00068 }
00069                 
00070 void CandOneToOneDeltaRMatcher::produce( Event& evt, const EventSetup& es ) {
00071   
00072   Handle<CandidateView> source;  
00073   Handle<CandidateView> matched;  
00074   evt.getByLabel( source_, source ) ;
00075   evt.getByLabel( matched_, matched ) ;
00076  
00077   edm::LogVerbatim("CandOneToOneDeltaRMatcher") << "======== Source Collection =======";
00078   for( CandidateView::const_iterator c = source->begin(); c != source->end(); ++c ) {
00079     edm::LogVerbatim("CandOneToOneDeltaRMatcher") << " pt source  " << c->pt() << " " << c->eta() << " " << c->phi()  << endl;
00080   }    
00081   edm::LogVerbatim("CandOneToOneDeltaRMatcher") << "======== Matched Collection =======";
00082   for( CandidateView::const_iterator c = matched->begin(); c != matched->end(); ++c ) {
00083     edm::LogVerbatim("CandOneToOneDeltaRMatcher") << " pt source  " << c->pt() << " " << c->eta() << " " << c->phi()  << endl;
00084   } 
00085 
00086   const int nSrc = source->size();
00087   const int nMtc = matched->size();
00088 
00089   const int nMin = min( source->size() , matched->size() );
00090   const int nMax = max( source->size() , matched->size() );
00091   if( nMin < 1 ) return;
00092 
00093   if( nSrc <= nMtc ) {
00094     for(CandidateView::const_iterator iSr  = source->begin();
00095         iSr != source->end();
00096         iSr++) {
00097       vector <float> tempAllDist;
00098       for(CandidateView::const_iterator iMt  = matched->begin();
00099           iMt != matched->end();
00100           iMt++) { 
00101         tempAllDist.push_back(DeltaR( iSr->p4() , iMt->p4() ) );
00102       }
00103       AllDist.push_back(tempAllDist);
00104       tempAllDist.clear();
00105     } 
00106   } else {
00107     for(CandidateView::const_iterator iMt  = matched->begin();
00108         iMt != matched->end();
00109         iMt++) {
00110       vector <float> tempAllDist;
00111       for(CandidateView::const_iterator iSr  = source->begin();
00112           iSr != source->end();
00113           iSr++) { 
00114         tempAllDist.push_back(DeltaR( iSr->p4() , iMt->p4() ) );
00115       }
00116       AllDist.push_back(tempAllDist);
00117       tempAllDist.clear();
00118     } 
00119   }
00120   
00121   /*
00122   edm::LogVerbatim("CandOneToOneDeltaRMatcher") << "======== The DeltaR Matrix =======";
00123   for(int m0=0; m0<nMin; m0++) {
00124     //    for(int m1=0; m1<nMax; m1++) {
00125       edm::LogVerbatim("CandOneToOneDeltaRMatcher") << setprecision(2) << fixed << (m1 AllDist[m0][m1] ;
00126     //}
00127     edm::LogVerbatim("CandOneToOneDeltaRMatcher") << "\n"; 
00128   }
00129   */
00130   
00131   // Loop size if Brute Force
00132   int nLoopToDo = (int) ( TMath::Factorial(nMax) / TMath::Factorial(nMax - nMin) );
00133   edm::LogVerbatim("CandOneToOneDeltaRMatcher") << "nLoop:" << nLoopToDo << endl;
00134   edm::LogVerbatim("CandOneToOneDeltaRMatcher") << "Choosen Algo is:" << algoMethod_ ;
00135   vector<int> bestCB;
00136 
00137   // Algo is Brute Force
00138   if( algoMethod_ == "BruteForce") {
00139 
00140     bestCB = AlgoBruteForce(nMin,nMax);
00141 
00142   // Algo is Switch Method
00143   } else if( algoMethod_ == "SwitchMode" ) {
00144 
00145     bestCB = AlgoSwitchMethod(nMin,nMax);
00146 
00147   // Algo is Brute Force if nLoop < 10000
00148   } else if( algoMethod_ == "MixMode" ) {
00149 
00150     if( nLoopToDo < 10000 ) {
00151       bestCB = AlgoBruteForce(nMin,nMax);
00152     } else { 
00153       bestCB = AlgoSwitchMethod(nMin,nMax);
00154     } 
00155 
00156   } else {
00157     throw cms::Exception("OneToOne Constructor") << "wrong matching method in ParameterSet";
00158   }
00159 
00160   for(int i1=0; i1<nMin; i1++) edm::LogVerbatim("CandOneToOneDeltaRMatcher") << "min: " << i1 << " " << bestCB[i1] << " " << AllDist[i1][bestCB[i1]];
00161 
00162 /*
00163   auto_ptr<CandViewMatchMap> matchMapSrMt( new CandViewMatchMap( CandViewMatchMap::ref_type( CandidateRefProd( source  ),
00164                                                                                              CandidateRefProd( matched )  ) ) );
00165   auto_ptr<CandViewMatchMap> matchMapMtSr( new CandViewMatchMap( CandViewMatchMap::ref_type( CandidateRefProd( matched ),
00166                                                                                              CandidateRefProd( source  )  ) ) );
00167 */
00168 
00169   auto_ptr<CandViewMatchMap> matchMapSrMt( new CandViewMatchMap() );
00170   auto_ptr<CandViewMatchMap> matchMapMtSr( new CandViewMatchMap() );
00171 
00172   for( int c = 0; c != nMin; c ++ ) {
00173     if( source->size() <= matched->size() ) {
00174       matchMapSrMt->insert( source ->refAt(c         ), matched->refAt(bestCB[c] ) );
00175       matchMapMtSr->insert( matched->refAt(bestCB[c] ), source ->refAt(c         ) );
00176     } else {
00177       matchMapSrMt->insert( source ->refAt(bestCB[c] ), matched->refAt(c         ) );
00178       matchMapMtSr->insert( matched->refAt(c         ), source ->refAt(bestCB[c] ) );
00179     }
00180   }
00181 
00182 /*
00183   for( int c = 0; c != nMin; c ++ ) {
00184     if( source->size() <= matched->size() ) { 
00185       matchMapSrMt->insert( CandidateRef( source,  c         ), CandidateRef( matched, bestCB[c] ) ); 
00186       matchMapMtSr->insert( CandidateRef( matched, bestCB[c] ), CandidateRef( source, c          ) );
00187     } else {
00188       matchMapSrMt->insert( CandidateRef( source,  bestCB[c] ), CandidateRef( matched, c         ) );
00189       matchMapMtSr->insert( CandidateRef( matched, c         ), CandidateRef( source,  bestCB[c] ) );
00190     }
00191   }
00192 */
00193   evt.put( matchMapSrMt, "src2mtc" );
00194   evt.put( matchMapMtSr, "mtc2src" );
00195 
00196   AllDist.clear();
00197 }
00198 
00199 
00200 double CandOneToOneDeltaRMatcher::lenght(const vector<int>& best) {
00201   double myLenght=0;
00202   int row=0;
00203   for(vector<int>::const_iterator it=best.begin(); it!=best.end(); it++ ) {             
00204     myLenght+=AllDist[row][*it];
00205     row++;
00206   }
00207   return myLenght;
00208 }
00209 
00210 // this is the Brute Force Algorithm
00211 // All the possible combination are checked
00212 // The best one is always found
00213 // Be carefull when you have high values for nMin and nMax --> the combinatorial could explode!
00214 // Sum(DeltaR) is minimized -->
00215 // 0.1 - 0.2 - 1.0 - 1.5 is lower than
00216 // 0.1 - 0.2 - 0.3 - 3.0 
00217 // Which one do you prefer? --> BruteForce select always the first
00218 
00219 vector<int> CandOneToOneDeltaRMatcher::AlgoBruteForce( int nMin, int nMax ) {
00220 
00221   vector<int> ca;
00222   vector<int> cb;
00223   vector<int> bestCB;
00224   float totalDeltaR=0;
00225   float BestTotalDeltaR=1000;
00226 
00227   for(int i1=0; i1<nMax; i1++) ca.push_back(i1);
00228   for(int i1=0; i1<nMin; i1++) cb.push_back(i1);
00229 
00230   do
00231     {
00232       //do your processing on the new combination here
00233       for(int cnt=0;cnt<TMath::Factorial(nMin); cnt++)
00234         {
00235           totalDeltaR = lenght(cb);
00236           if ( totalDeltaR < BestTotalDeltaR ) {
00237             BestTotalDeltaR = totalDeltaR;
00238             bestCB=cb;
00239           }
00240           next_permutation( cb.begin() , cb.end() );
00241         }
00242     }
00243   while(next_combination( ca.begin() , ca.end() , cb.begin() , cb.end() ));
00244   
00245   return bestCB;
00246 }
00247 
00248 // This method (Developed originally by Daniele Benedetti) check for the best combination
00249 // choosing the minimum DeltaR for each line in AllDist matrix
00250 // If no repeated row is found: ie (line,col)=(1,3) and (2,3) --> same as BruteForce
00251 // If repetition --> set the higher DeltaR between  the 2 repetition to 1000 and re-check best combination
00252 // Iterate until no repetition  
00253 // No guaranted minimum for Sum(DeltaR)
00254 // If you have:
00255 // 0.1 - 0.2 - 1.0 - 1.5 is lower than
00256 // 0.1 - 0.2 - 0.3 - 3.0 
00257 // SwitchMethod normally select the second solution
00258 
00259 vector<int> CandOneToOneDeltaRMatcher::AlgoSwitchMethod( int nMin, int nMax ) {
00260 
00261   vector<int> bestCB;
00262   for(int i1=0; i1<nMin; i1++) {
00263     int minInd=0;
00264     for(int i2=1; i2<nMax; i2++) if( AllDist[i1][i2] < AllDist[i1][minInd] ) minInd = i2; 
00265     bestCB.push_back(minInd);
00266   }
00267 
00268   bool inside = true;
00269   while( inside ) {
00270     inside = false;
00271     for(int i1=0;i1<nMin;i1++){
00272       for(int i2=i1+1;i2<nMin;i2++){
00273         if ( bestCB[i1] == bestCB[i2] ) {
00274           inside = true;
00275           if ( AllDist[i1][(bestCB[i1])] <= AllDist[i2][(bestCB[i2])]) {
00276             AllDist[i2][(bestCB[i2])]= 1000;
00277             int minInd=0;
00278             for(int i3=1; i3<nMax; i3++) if( AllDist[i2][i3] < AllDist[i2][minInd] ) minInd = i3; 
00279             bestCB[i2]= minInd;
00280           }  else {
00281             AllDist[i1][(bestCB[i1])]= 1000;
00282             int minInd=0;
00283             for(int i3=1; i3<nMax; i3++) if( AllDist[i1][i3] < AllDist[i1][minInd] ) minInd = i3; 
00284             bestCB[i1]= minInd;
00285           }
00286         } // End if
00287       } 
00288     }
00289   } // End while
00290 
00291   return bestCB;
00292 
00293 }
00294 
00295 #include "FWCore/PluginManager/interface/ModuleDef.h"
00296 #include "FWCore/Framework/interface/MakerMacros.h"
00297 
00298 DEFINE_FWK_MODULE( CandOneToOneDeltaRMatcher );
/afs/cern.ch/work/a/aaltunda/public/www/CMSSW_6_2_7/src/PhysicsTools/JetMCAlgos/plugins/CandOneToOneDeltaRMatcher.cc
