CandOneToOneDeltaRMatcher.cc

Go to the documentation of this file.

/* \class CandOneToOneDeltaRMatcher
 *
 * Producer for simple match map
 * to match two collections of candidate
 * with one-to-One matching
 * minimizing Sum(DeltaR)
 *
 */

#include "FWCore/Framework/interface/EDProducer.h"
#include "FWCore/ParameterSet/interface/ParameterSetfwd.h"

#include "DataFormats/Candidate/interface/Candidate.h"

#include<vector>
#include<iostream>

class CandOneToOneDeltaRMatcher : public edm::EDProducer {
 public:
  CandOneToOneDeltaRMatcher( const edm::ParameterSet & );
  ~CandOneToOneDeltaRMatcher() override;
 private:
  void produce( edm::Event&, const edm::EventSetup& ) override;
  double lenght( const std::vector<int>& );
  std::vector<int> AlgoBruteForce(int, int);
  std::vector<int> AlgoSwitchMethod(int, int);

  edm::EDGetTokenT<reco::CandidateView> sourceToken_;
  edm::EDGetTokenT<reco::CandidateView> matchedToken_;
  std::vector < std::vector<float> > AllDist;
  std::string algoMethod_;

};

#include "PhysicsTools/JetMCUtils/interface/combination.h"

#include "FWCore/Framework/interface/ESHandle.h"
#include "FWCore/Framework/interface/Event.h"
#include "FWCore/Framework/interface/EventSetup.h"
#include "FWCore/Utilities/interface/InputTag.h"
#include "FWCore/Utilities/interface/EDMException.h"
#include "FWCore/ParameterSet/interface/ParameterSet.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"

#include "DataFormats/Common/interface/Handle.h"
#include "DataFormats/Candidate/interface/LeafCandidate.h"
#include "DataFormats/Candidate/interface/CandMatchMap.h"
#include "DataFormats/Candidate/interface/CandidateFwd.h"


#include <Math/VectorUtil.h>
#include <TMath.h>

using namespace edm;
using namespace std;
using namespace reco;
using namespace ROOT::Math::VectorUtil;
using namespace stdcomb;

CandOneToOneDeltaRMatcher::CandOneToOneDeltaRMatcher( const ParameterSet & cfg ) :
  sourceToken_( consumes<CandidateView>( cfg.getParameter<InputTag>( "src" ) ) ),
  matchedToken_( consumes<CandidateView>( cfg.getParameter<InputTag>( "matched" ) ) ),
  algoMethod_( cfg.getParameter<string>( "algoMethod" ) ) {
  produces<CandViewMatchMap>("src2mtc");
  produces<CandViewMatchMap>("mtc2src");
}

CandOneToOneDeltaRMatcher::~CandOneToOneDeltaRMatcher() {
}

void CandOneToOneDeltaRMatcher::produce( Event& evt, const EventSetup& es ) {

  Handle<CandidateView> source;
  Handle<CandidateView> matched;
  evt.getByToken( sourceToken_, source ) ;
  evt.getByToken( matchedToken_, matched ) ;

  edm::LogVerbatim("CandOneToOneDeltaRMatcher") << "======== Source Collection =======";
  for( CandidateView::const_iterator c = source->begin(); c != source->end(); ++c ) {
    edm::LogVerbatim("CandOneToOneDeltaRMatcher") << " pt source  " << c->pt() << " " << c->eta() << " " << c->phi()  << endl;
  }
  edm::LogVerbatim("CandOneToOneDeltaRMatcher") << "======== Matched Collection =======";
  for( CandidateView::const_iterator c = matched->begin(); c != matched->end(); ++c ) {
    edm::LogVerbatim("CandOneToOneDeltaRMatcher") << " pt source  " << c->pt() << " " << c->eta() << " " << c->phi()  << endl;
  }

  const int nSrc = source->size();
  const int nMtc = matched->size();

  const int nMin = min( source->size() , matched->size() );
  const int nMax = max( source->size() , matched->size() );
  if( nMin < 1 ) return;

  if( nSrc <= nMtc ) {
    for(CandidateView::const_iterator iSr  = source->begin();
    iSr != source->end();
    iSr++) {
      vector <float> tempAllDist;
      for(CandidateView::const_iterator iMt  = matched->begin();
      iMt != matched->end();
      iMt++) {
    tempAllDist.push_back(DeltaR( iSr->p4() , iMt->p4() ) );
      }
      AllDist.push_back(tempAllDist);
      tempAllDist.clear();
    }
  } else {
    for(CandidateView::const_iterator iMt  = matched->begin();
    iMt != matched->end();
    iMt++) {
      vector <float> tempAllDist;
      for(CandidateView::const_iterator iSr  = source->begin();
      iSr != source->end();
      iSr++) {
    tempAllDist.push_back(DeltaR( iSr->p4() , iMt->p4() ) );
      }
      AllDist.push_back(tempAllDist);
      tempAllDist.clear();
    }
  }

  /*
  edm::LogVerbatim("CandOneToOneDeltaRMatcher") << "======== The DeltaR Matrix =======";
  for(int m0=0; m0<nMin; m0++) {
    //    for(int m1=0; m1<nMax; m1++) {
      edm::LogVerbatim("CandOneToOneDeltaRMatcher") << setprecision(2) << fixed << (m1 AllDist[m0][m1] ;
    //}
    edm::LogVerbatim("CandOneToOneDeltaRMatcher") << "\n";
  }
  */

  // Loop size if Brute Force
  int nLoopToDo = (int) ( TMath::Factorial(nMax) / TMath::Factorial(nMax - nMin) );
  edm::LogVerbatim("CandOneToOneDeltaRMatcher") << "nLoop:" << nLoopToDo << endl;
  edm::LogVerbatim("CandOneToOneDeltaRMatcher") << "Choosen Algo is:" << algoMethod_ ;
  vector<int> bestCB;

  // Algo is Brute Force
  if( algoMethod_ == "BruteForce") {

    bestCB = AlgoBruteForce(nMin,nMax);

  // Algo is Switch Method
  } else if( algoMethod_ == "SwitchMode" ) {

    bestCB = AlgoSwitchMethod(nMin,nMax);

  // Algo is Brute Force if nLoop < 10000
  } else if( algoMethod_ == "MixMode" ) {

    if( nLoopToDo < 10000 ) {
      bestCB = AlgoBruteForce(nMin,nMax);
    } else {
      bestCB = AlgoSwitchMethod(nMin,nMax);
    }

  } else {
    throw cms::Exception("OneToOne Constructor") << "wrong matching method in ParameterSet";
  }

  for(int i1=0; i1<nMin; i1++) edm::LogVerbatim("CandOneToOneDeltaRMatcher") << "min: " << i1 << " " << bestCB[i1] << " " << AllDist[i1][bestCB[i1]];

/*
  auto matchMapSrMt = std::make_unique<CandViewMatchMap>(CandViewMatchMap::ref_type( CandidateRefProd( source  ),
                                                                                             CandidateRefProd( matched ) ) );
  auto matchMapMtSr = std::make_unique<CandViewMatchMap>(CandViewMatchMap::ref_type( CandidateRefProd( matched ),
                                                                                             CandidateRefProd( source ) ) );
*/

  auto matchMapSrMt = std::make_unique<CandViewMatchMap>();
  auto matchMapMtSr = std::make_unique<CandViewMatchMap>();

  for( int c = 0; c != nMin; c ++ ) {
    if( source->size() <= matched->size() ) {
      matchMapSrMt->insert( source ->refAt(c         ), matched->refAt(bestCB[c] ) );
      matchMapMtSr->insert( matched->refAt(bestCB[c] ), source ->refAt(c         ) );
    } else {
      matchMapSrMt->insert( source ->refAt(bestCB[c] ), matched->refAt(c         ) );
      matchMapMtSr->insert( matched->refAt(c         ), source ->refAt(bestCB[c] ) );
    }
  }

/*
  for( int c = 0; c != nMin; c ++ ) {
    if( source->size() <= matched->size() ) {
      matchMapSrMt->insert( CandidateRef( source,  c         ), CandidateRef( matched, bestCB[c] ) );
      matchMapMtSr->insert( CandidateRef( matched, bestCB[c] ), CandidateRef( source, c          ) );
    } else {
      matchMapSrMt->insert( CandidateRef( source,  bestCB[c] ), CandidateRef( matched, c         ) );
      matchMapMtSr->insert( CandidateRef( matched, c         ), CandidateRef( source,  bestCB[c] ) );
    }
  }
*/
  evt.put(std::move(matchMapSrMt), "src2mtc" );
  evt.put(std::move(matchMapMtSr), "mtc2src" );

  AllDist.clear();
}


double CandOneToOneDeltaRMatcher::lenght(const vector<int>& best) {
  double myLenght=0;
  int row=0;
  for(vector<int>::const_iterator it=best.begin(); it!=best.end(); it++ ) {
    myLenght+=AllDist[row][*it];
    row++;
  }
  return myLenght;
}

// this is the Brute Force Algorithm
// All the possible combination are checked
// The best one is always found
// Be carefull when you have high values for nMin and nMax --> the combinatorial could explode!
// Sum(DeltaR) is minimized -->
// 0.1 - 0.2 - 1.0 - 1.5 is lower than
// 0.1 - 0.2 - 0.3 - 3.0
// Which one do you prefer? --> BruteForce select always the first

vector<int> CandOneToOneDeltaRMatcher::AlgoBruteForce( int nMin, int nMax ) {

  vector<int> ca;
  vector<int> cb;
  vector<int> bestCB;
  float totalDeltaR=0;
  float BestTotalDeltaR=1000;

  for(int i1=0; i1<nMax; i1++) ca.push_back(i1);
  for(int i1=0; i1<nMin; i1++) cb.push_back(i1);

  do
    {
      //do your processing on the new combination here
      for(int cnt=0;cnt<TMath::Factorial(nMin); cnt++)
    {
      totalDeltaR = lenght(cb);
      if ( totalDeltaR < BestTotalDeltaR ) {
        BestTotalDeltaR = totalDeltaR;
        bestCB=cb;
      }
      next_permutation( cb.begin() , cb.end() );
    }
    }
  while(next_combination( ca.begin() , ca.end() , cb.begin() , cb.end() ));

  return bestCB;
}

// This method (Developed originally by Daniele Benedetti) check for the best combination
// choosing the minimum DeltaR for each line in AllDist matrix
// If no repeated row is found: ie (line,col)=(1,3) and (2,3) --> same as BruteForce
// If repetition --> set the higher DeltaR between  the 2 repetition to 1000 and re-check best combination
// Iterate until no repetition
// No guaranted minimum for Sum(DeltaR)
// If you have:
// 0.1 - 0.2 - 1.0 - 1.5 is lower than
// 0.1 - 0.2 - 0.3 - 3.0
// SwitchMethod normally select the second solution

vector<int> CandOneToOneDeltaRMatcher::AlgoSwitchMethod( int nMin, int nMax ) {

  vector<int> bestCB;
  for(int i1=0; i1<nMin; i1++) {
    int minInd=0;
    for(int i2=1; i2<nMax; i2++) if( AllDist[i1][i2] < AllDist[i1][minInd] ) minInd = i2;
    bestCB.push_back(minInd);
  }

  bool inside = true;
  while( inside ) {
    inside = false;
    for(int i1=0;i1<nMin;i1++){
      for(int i2=i1+1;i2<nMin;i2++){
    if ( bestCB[i1] == bestCB[i2] ) {
      inside = true;
      if ( AllDist[i1][(bestCB[i1])] <= AllDist[i2][(bestCB[i2])]) {
        AllDist[i2][(bestCB[i2])]= 1000;
        int minInd=0;
        for(int i3=1; i3<nMax; i3++) if( AllDist[i2][i3] < AllDist[i2][minInd] ) minInd = i3;
        bestCB[i2]= minInd;
      }  else {
        AllDist[i1][(bestCB[i1])]= 1000;
        int minInd=0;
        for(int i3=1; i3<nMax; i3++) if( AllDist[i1][i3] < AllDist[i1][minInd] ) minInd = i3;
        bestCB[i1]= minInd;
      }
    } // End if
      }
    }
  } // End while

  return bestCB;

}

#include "FWCore/PluginManager/interface/ModuleDef.h"
#include "FWCore/Framework/interface/MakerMacros.h"

DEFINE_FWK_MODULE( CandOneToOneDeltaRMatcher );