#include <JetPartonMatching.h>

Public Types
enum	algorithms { totalMinDist, minSumDist, ptOrderedMinDist, unambiguousOnly }
typedef std::vector< std::pair < unsigned int, int > >	MatchingCollection
Public Member Functions
double	getDistanceForParton (const unsigned int part, const unsigned int comb=0)
std::vector< int >	getMatchesForPartons (const unsigned int comb=0)
int	getMatchForParton (const unsigned int part, const unsigned int comb=0)
unsigned int	getNumberOfAvailableCombinations ()
int	getNumberOfUnmatchedPartons (const unsigned int comb=0)
double	getSumDeltaE (const unsigned int comb=0)
double	getSumDeltaPt (const unsigned int comb=0)
double	getSumDeltaR (const unsigned int comb=0)
double	getSumDistances (const unsigned int comb=0)
	JetPartonMatching (const std::vector< const reco::Candidate * > &, const std::vector< pat::Jet > &, const int, const bool, const bool, const double)
	JetPartonMatching (const std::vector< const reco::Candidate * > &, const std::vector< const reco::Candidate * > &, const int, const bool, const bool, const double)
	JetPartonMatching (const std::vector< const reco::Candidate * > &, const std::vector< reco::CaloJet > &, const int, const bool, const bool, const double)
	JetPartonMatching ()
	JetPartonMatching (const std::vector< const reco::Candidate * > &, const std::vector< reco::GenJet > &, const int, const bool, const bool, const double)
void	print ()
	~JetPartonMatching ()
Private Member Functions
void	calculate ()
double	distance (const math::XYZTLorentzVector &, const math::XYZTLorentzVector &)
void	matchingMinSumDist ()
void	matchingPtOrderedMinDist ()
void	matchingTotalMinDist ()
void	matchingUnambiguousOnly ()
void	minSumDist_recursion (const unsigned int, std::vector< unsigned int > &, std::vector< bool > &, std::vector< std::pair< double, MatchingCollection > > &)
Private Attributes
int	algorithm_
std::vector< const reco::Candidate * >	jets
std::vector< MatchingCollection >	matching
double	maxDist_
std::vector< unsigned int >	numberOfUnmatchedPartons
std::vector< const reco::Candidate * >	partons
std::vector< double >	sumDeltaE
std::vector< double >	sumDeltaPt
std::vector< double >	sumDeltaR
bool	useDeltaR_
bool	useMaxDist_

Detailed Description

Definition at line 11 of file JetPartonMatching.h.

Member Typedef Documentation

typedef std::vector< std::pair<unsigned int, int> > JetPartonMatching::MatchingCollection

Definition at line 16 of file JetPartonMatching.h.

Member Enumeration Documentation

enum JetPartonMatching::algorithms

Enumerator:

totalMinDist
minSumDist
ptOrderedMinDist
unambiguousOnly

Definition at line 17 of file JetPartonMatching.h.

{ totalMinDist, minSumDist, ptOrderedMinDist, unambiguousOnly };

Constructor & Destructor Documentation

JetPartonMatching::JetPartonMatching ( ) [inline]

Definition at line 19 of file JetPartonMatching.h.

{};

JetPartonMatching::JetPartonMatching	(	const std::vector< const reco::Candidate * > &	p,
		const std::vector< reco::GenJet > &	j,
		const int	algorithm = `totalMinDist`,
		const bool	useMaxDist = `true`,
		const bool	useDeltaR = `true`,
		const double	maxDist = `0.3`
	)

Definition at line 5 of file JetPartonMatching.cc.

References calculate(), i, and jets.

  : partons(p), algorithm_(algorithm), useMaxDist_(useMaxDist), useDeltaR_(useDeltaR), maxDist_(maxDist)
{
  std::vector<const reco::Candidate*> js;
  for(unsigned int i=0; i<j.size(); ++i) js.push_back( &(j[i]) );
  jets = js;
  calculate();
}

JetPartonMatching::JetPartonMatching	(	const std::vector< const reco::Candidate * > &	p,
		const std::vector< reco::CaloJet > &	j,
		const int	algorithm = `totalMinDist`,
		const bool	useMaxDist = `true`,
		const bool	useDeltaR = `true`,
		const double	maxDist = `0.3`
	)

Definition at line 16 of file JetPartonMatching.cc.

References calculate(), i, and jets.

  : partons(p), algorithm_(algorithm), useMaxDist_(useMaxDist), useDeltaR_(useDeltaR), maxDist_(maxDist)
{
  std::vector<const reco::Candidate*> js;
  for(unsigned int i = 0; i < j.size(); ++i) js.push_back( &(j[i]) );
  jets = js; 
  calculate(); 
}

JetPartonMatching::JetPartonMatching	(	const std::vector< const reco::Candidate * > &	p,
		const std::vector< pat::Jet > &	j,
		const int	algorithm = `totalMinDist`,
		const bool	useMaxDist = `true`,
		const bool	useDeltaR = `true`,
		const double	maxDist = `0.3`
	)

Definition at line 27 of file JetPartonMatching.cc.

References calculate(), i, and jets.

  : partons(p), algorithm_(algorithm), useMaxDist_(useMaxDist), useDeltaR_(useDeltaR), maxDist_(maxDist)
{
  std::vector<const reco::Candidate*> js;
  for(unsigned int i = 0; i < j.size(); ++i) js.push_back( &(j[i]) );
  jets = js; 
  calculate(); 
}

JetPartonMatching::JetPartonMatching	(	const std::vector< const reco::Candidate * > &	p,
		const std::vector< const reco::Candidate * > &	j,
		const int	algorithm = `totalMinDist`,
		const bool	useMaxDist = `true`,
		const bool	useDeltaR = `true`,
		const double	maxDist = `0.3`
	)

Definition at line 38 of file JetPartonMatching.cc.

References calculate().

  : partons(p), jets(j), algorithm_(algorithm), useMaxDist_(useMaxDist), useDeltaR_(useDeltaR), maxDist_(maxDist)
{
  calculate();
}

JetPartonMatching::~JetPartonMatching ( ) [inline]

Definition at line 28 of file JetPartonMatching.h.

{};

Member Function Documentation

void JetPartonMatching::calculate ( ) [private]

Definition at line 47 of file JetPartonMatching.cc.

References algorithm_, bookConverter::comb, distance(), relval_parameters_module::energy, first, getMatchForParton(), i, jets, match(), matching, matchingMinSumDist(), matchingPtOrderedMinDist(), matchingTotalMinDist(), matchingUnambiguousOnly(), minSumDist, numberOfUnmatchedPartons, p4, partons, benchmark_cfg::pdgId, ptOrderedMinDist, edm::second(), python::multivaluedict::sort(), sumDeltaE, sumDeltaPt, sumDeltaR, totalMinDist, unambiguousOnly, and useMaxDist_.

Referenced by JetPartonMatching().

{
  // use maximal distance between objects 
  // in case of unambiguousOnly algorithmm
  if(algorithm_==unambiguousOnly) useMaxDist_=true;

  // check if there are empty partons in
  // the vector, which happpens if the 
  // event is not ttbar or the decay is 
  // not as expected
  bool emptyParton=false;
  for(unsigned int ip=0; ip<partons.size(); ++ip){
    if( partons[ip]->pdgId() ==0 ){
      emptyParton=true;
      break;
    }
  }

  // switch algorithm, default is to match
  // on the minimal sum of the distance 
  // (if jets or a parton is empty fill match with blanks)
  if( jets.empty() || emptyParton ) {
    MatchingCollection dummyMatch;
    for(unsigned int ip=0; ip<partons.size(); ++ip)
      dummyMatch.push_back( std::make_pair(ip, -1) );
    matching.push_back( dummyMatch );
  }
  else {
    switch(algorithm_) {
      
    case totalMinDist: 
      matchingTotalMinDist();    
      break;
      
    case minSumDist: 
      matchingMinSumDist();
      break;
      
    case ptOrderedMinDist: 
      matchingPtOrderedMinDist();
      break;
      
    case unambiguousOnly:
      matchingUnambiguousOnly();
      break;
      
    default:
      matchingMinSumDist();
    }
  }

  numberOfUnmatchedPartons.clear();
  sumDeltaE .clear();
  sumDeltaPt.clear();
  sumDeltaR .clear();
  for(unsigned int comb = 0; comb < matching.size(); ++comb) {
    MatchingCollection match = matching[comb];
    std::sort(match.begin(), match.end());
    matching[comb] = match;
    
    int nUnmatchedPartons = partons.size();
    for(unsigned int part=0; part<partons.size(); ++part)
      if(getMatchForParton(part,comb)>=0) --nUnmatchedPartons;

    double sumDE  = -999.;
    double sumDPt = -999.;
    double sumDR  = -999.;
    if(nUnmatchedPartons==0){
      sumDE  = 0;
      sumDPt = 0;
      sumDR  = 0;
      for(unsigned int i=0; i<match.size(); ++i){
        sumDE  += fabs(partons[match[i].first]->energy() - jets[match[i].second]->energy());
        sumDPt += fabs(partons[match[i].first]->pt()     - jets[match[i].second]->pt());
        sumDR  += distance(partons[match[i].first]->p4(), jets[match[i].second]->p4());
      }
    }

    numberOfUnmatchedPartons.push_back( nUnmatchedPartons );
    sumDeltaE .push_back( sumDE  );
    sumDeltaPt.push_back( sumDPt );
    sumDeltaR .push_back( sumDR  );
  }
}

double JetPartonMatching::distance	(	const math::XYZTLorentzVector &	v1,
		const math::XYZTLorentzVector &	v2
	)		`[private]`

Definition at line 133 of file JetPartonMatching.cc.

References useDeltaR_.

Referenced by calculate(), getDistanceForParton(), matchingPtOrderedMinDist(), matchingTotalMinDist(), matchingUnambiguousOnly(), and minSumDist_recursion().

{
  // calculate the distance between two lorentz vectors 
  // using DeltaR(eta, phi) or normal space angle(theta, phi)
  if(useDeltaR_) return ROOT::Math::VectorUtil::DeltaR(v1, v2);
  return ROOT::Math::VectorUtil::Angle(v1, v2);
}

double JetPartonMatching::getDistanceForParton	(	const unsigned int	part,
		const unsigned int	comb = `0`
	)

Definition at line 358 of file JetPartonMatching.cc.

References distance(), getMatchForParton(), jets, p4, and partons.

Referenced by getSumDistances(), TtSemiEvtSolutionMaker::produce(), and TtHadEvtSolutionMaker::produce().

{
  // get the distance between parton and its best matched jet
  if(getMatchForParton(part, comb) < 0) return -999.;
  return distance( jets[getMatchForParton(part,comb)]->p4(), partons[part]->p4() );
}

std::vector< int > JetPartonMatching::getMatchesForPartons ( const unsigned int comb = 0 )

Definition at line 347 of file JetPartonMatching.cc.

References getMatchForParton(), and partons.

Referenced by TtJetPartonMatch< C >::produce().

{
  // return a vector with the indices of the matched jets
  // (ordered according to the vector of partons)
  std::vector<int> jetIndices;
  for(unsigned int part=0; part<partons.size(); ++part)
    jetIndices.push_back( getMatchForParton(part, comb) );
  return jetIndices;
}

int JetPartonMatching::getMatchForParton	(	const unsigned int	part,
		const unsigned int	comb = `0`
	)

Definition at line 336 of file JetPartonMatching.cc.

References matching, and findQualityFiles::size.

Referenced by calculate(), getDistanceForParton(), getMatchesForPartons(), print(), TtSemiEvtSolutionMaker::produce(), and TtHadEvtSolutionMaker::produce().

{
  // return index of the matched jet for a given parton
  // (if arguments for parton index and combinatoric index
  // are in the valid range)
  if(comb >= matching.size()) return -9;
  if(part >= matching[comb].size()) return -9;
  return (matching[comb])[part].second;
}

unsigned int JetPartonMatching::getNumberOfAvailableCombinations ( ) [inline]

Definition at line 31 of file JetPartonMatching.h.

References matching.

Referenced by print(), and TtJetPartonMatch< C >::produce().

{ return matching.size(); }

int JetPartonMatching::getNumberOfUnmatchedPartons ( const unsigned int comb = 0 ) [inline]

Definition at line 32 of file JetPartonMatching.h.

References bookConverter::comb, and numberOfUnmatchedPartons.

{ return (comb<numberOfUnmatchedPartons.size() ? (int)numberOfUnmatchedPartons[comb] : -1 ); }

double JetPartonMatching::getSumDeltaE ( const unsigned int comb = 0 ) [inline]

Definition at line 39 of file JetPartonMatching.h.

References bookConverter::comb, findQualityFiles::size, and sumDeltaE.

Referenced by print().

{ return (comb<sumDeltaE .size() ? sumDeltaE [comb] : -999.); }

double JetPartonMatching::getSumDeltaPt ( const unsigned int comb = 0 ) [inline]

Definition at line 40 of file JetPartonMatching.h.

References bookConverter::comb, and sumDeltaPt.

Referenced by print(), and TtJetPartonMatch< C >::produce().

{ return (comb<sumDeltaPt.size() ? sumDeltaPt[comb] : -999.); }

double JetPartonMatching::getSumDeltaR ( const unsigned int comb = 0 ) [inline]

Definition at line 41 of file JetPartonMatching.h.

References bookConverter::comb, findQualityFiles::size, and sumDeltaR.

Referenced by print(), and TtJetPartonMatch< C >::produce().

{ return (comb<sumDeltaR .size() ? sumDeltaR [comb] : -999.); }

double JetPartonMatching::getSumDistances ( const unsigned int comb = 0 )

Definition at line 366 of file JetPartonMatching.cc.

References getDistanceForParton(), and partons.

Referenced by TtSemiEvtSolutionMaker::produce(), and TtHadEvtSolutionMaker::produce().

{
  // get sum of distances between partons and matched jets
  double sumDists = 0.;
  for(unsigned int part=0; part<partons.size(); ++part){
    double dist = getDistanceForParton(part, comb);
    if(dist < 0.) return -999.;
    sumDists += dist;
  }
  return sumDists;
}

void JetPartonMatching::matchingMinSumDist ( ) [private]

Definition at line 221 of file JetPartonMatching.cc.

References i, jets, matching, minSumDist_recursion(), partons, and python::multivaluedict::sort().

Referenced by calculate().

{
  // match partons to jets with minimal sum of
  // the distances between all partons and jets

  std::vector<std::pair<double, MatchingCollection> > distMatchVec;

  std::vector<bool> usedJets;
  for(unsigned int i=0; i<jets.size(); ++i){
    usedJets.push_back(false);
  }

  std::vector<unsigned int> jetIndices;
  jetIndices.reserve(partons.size());

  minSumDist_recursion(0, jetIndices, usedJets, distMatchVec);

  std::sort(distMatchVec.begin(), distMatchVec.end());

  matching.clear();

  if(distMatchVec.empty()) {
    MatchingCollection dummyMatch;
    for(unsigned int ip=0; ip<partons.size(); ++ip)
      dummyMatch.push_back(std::make_pair(ip, -1));
    matching.push_back( dummyMatch );
  }
  else
    for(unsigned int i=0; i<distMatchVec.size(); ++i)
      matching.push_back( distMatchVec[i].second );

  return;
}

void JetPartonMatching::matchingPtOrderedMinDist ( ) [private]

Definition at line 256 of file JetPartonMatching.cc.

References distance(), jets, match(), matching, maxDist_, p4, partons, edm::second(), python::multivaluedict::sort(), and useMaxDist_.

Referenced by calculate().

{
  // match partons to jets with minimal sum of
  // the distances between all partons and jets
  // order partons in pt first
  std::vector<std::pair <double, unsigned int> > ptOrderedPartons;

  for(unsigned int ip=0; ip<partons.size(); ++ip)
    ptOrderedPartons.push_back(std::make_pair(partons[ip]->pt(), ip));

  std::sort(ptOrderedPartons.begin(), ptOrderedPartons.end());
  std::reverse(ptOrderedPartons.begin(), ptOrderedPartons.end());

  std::vector<unsigned int> jetIndices;
  for(unsigned int ij=0; ij<jets.size(); ++ij) jetIndices.push_back(ij);

  MatchingCollection match;

  for(unsigned int ip=0; ip<ptOrderedPartons.size(); ++ip){
    double minDist = 999.;
    int ijMin = -1;

    for(unsigned int ij=0; ij<jetIndices.size(); ++ij){
      double dist = distance(partons[ptOrderedPartons[ip].second]->p4(), jets[jetIndices[ij]]->p4());
      if(dist < minDist){
        if(!useMaxDist_ || dist <= maxDist_){
          minDist = dist;
          ijMin = ij;
        }
      }
    }
    
    if(ijMin >= 0){
      match.push_back( std::make_pair(ptOrderedPartons[ip].second, jetIndices[ijMin]) );
      jetIndices.erase(jetIndices.begin() + ijMin, jetIndices.begin() + ijMin + 1);
    }
    else
      match.push_back( std::make_pair(ptOrderedPartons[ip].second, -1) );
  }

  matching.clear();
  matching.push_back( match );
  return;
}

void JetPartonMatching::matchingTotalMinDist ( ) [private]

Definition at line 142 of file JetPartonMatching.cc.

References a, distance(), first, jets, match(), matching, maxDist_, p4, partons, python::multivaluedict::sort(), and useMaxDist_.

Referenced by calculate().

{
  // match parton to jet with shortest distance
  // starting with the shortest distance available
  // apply some outlier rejection if desired

  // prepare vector of pairs with distances between
  // all partons to all jets in the input vectors
  std::vector< std::pair<double, unsigned int> > distances;
  for(unsigned int ip=0; ip<partons.size(); ++ip){
    for(unsigned int ij=0; ij<jets.size(); ++ij){ 
      double dist = distance(jets[ij]->p4(), partons[ip]->p4());
      distances.push_back(std::pair<double, unsigned int>(dist, ip*jets.size()+ij));
    }
  }
  std::sort(distances.begin(), distances.end());

  MatchingCollection match;

  while(match.size() < partons.size()){
    unsigned int partonIndex = distances[0].second/jets.size();
    int jetIndex = distances[0].second-jets.size()*partonIndex;
    
    // use primitive outlier rejection if desired
    if(useMaxDist_&& distances[0].first>maxDist_) jetIndex = -1;

    // prevent underflow in case of too few jets
    if( distances.empty() )
      match.push_back(std::make_pair(partonIndex, -1));
    else
      match.push_back(std::make_pair(partonIndex, jetIndex));
    
    // remove all values for the matched parton 
    // and the matched jet
    for(unsigned int a=0; a<distances.size(); ++a){
      unsigned int pIndex = distances[a].second/jets.size();
      int jIndex = distances[a].second-jets.size()*pIndex;
      if((pIndex == partonIndex) || (jIndex == jetIndex)){
        distances.erase(distances.begin()+a, distances.begin()+a+1); 
        --a;
      }
    }
  }

  matching.clear();
  matching.push_back( match );
  return;
}

void JetPartonMatching::matchingUnambiguousOnly ( ) [private]

Definition at line 302 of file JetPartonMatching.cc.

References distance(), jets, match(), matching, maxDist_, p4, and partons.

Referenced by calculate().

{
  // match partons to jets, only accept event 
  // if there are no ambiguities
  std::vector<bool> jetMatched;
  for(unsigned int ij=0; ij<jets.size(); ++ij) jetMatched.push_back(false);

  MatchingCollection match;
  
  for(unsigned int ip=0; ip<partons.size(); ++ip){
    int iMatch = -1;
    for(unsigned int ij=0; ij<jets.size(); ++ij){
      double dist = distance(partons[ip]->p4(), jets[ij]->p4());
      if(dist <= maxDist_){
        if(!jetMatched[ij]){ // new match for jet
          jetMatched[ij] = true;
          if(iMatch == -1) // new match for parton and jet
            iMatch = ij;
          else // parton already matched: ambiguity!
            iMatch = -2;
        }
        else // jet already matched: ambiguity!
          iMatch = -2;
      }
    }
    match.push_back(std::make_pair(ip, iMatch));
  }

  matching.clear();
  matching.push_back( match );
  return;
}

void JetPartonMatching::minSumDist_recursion	(	const unsigned int	ip,
		std::vector< unsigned int > &	jetIndices,
		std::vector< bool > &	usedJets,
		std::vector< std::pair< double, MatchingCollection > > &	distMatchVec
	)		`[private]`

Definition at line 192 of file JetPartonMatching.cc.

References distance(), jets, match(), maxDist_, p4, partons, and useMaxDist_.

Referenced by matchingMinSumDist().

{
  // build up jet combinations recursively
  if(ip<partons.size()){
    for(unsigned int ij=0; ij<jets.size(); ++ij){
      if(usedJets[ij]) continue;
      usedJets[ij] = true;
      jetIndices[ip] = ij;
      minSumDist_recursion(ip+1, jetIndices, usedJets, distMatchVec);
      usedJets[ij] = false;
    }
    return;
  }

  // calculate sumDist for each completed combination
  double sumDist = 0;
  MatchingCollection match;
  for(unsigned int ip=0; ip<partons.size(); ++ip){
    double dist  = distance(partons[ip]->p4(), jets[jetIndices[ip]]->p4());
    if(useMaxDist_ && dist > maxDist_) return; // outlier rejection
    sumDist += distance(partons[ip]->p4(), jets[jetIndices[ip]]->p4());
    match.push_back(std::make_pair(ip, jetIndices[ip]));
  }

  distMatchVec.push_back( std::make_pair(sumDist, match)  );
  return;
}

void JetPartonMatching::print ( void )

Definition at line 379 of file JetPartonMatching.cc.

References algorithm_, bookConverter::comb, getMatchForParton(), getNumberOfAvailableCombinations(), getSumDeltaE(), getSumDeltaPt(), getSumDeltaR(), jets, funct::log(), matching, maxDist_, minSumDist, partons, ptOrderedMinDist, totalMinDist, unambiguousOnly, useDeltaR_, and useMaxDist_.

Referenced by TtJetPartonMatch< C >::produce().

{
  //report using MessageLogger
  edm::LogInfo log("JetPartonMatching");
  log << "++++++++++++++++++++++++++++++++++++++++++++++ \n";
  log << " algorithm : ";
  switch(algorithm_) {
  case totalMinDist     : log << "totalMinDist    "; break;
  case minSumDist       : log << "minSumDist      "; break;
  case ptOrderedMinDist : log << "ptOrderedMinDist"; break;
  case unambiguousOnly  : log << "unambiguousOnly "; break;
  default               : log << "UNKNOWN         ";
  }
  log << "\n";
  log << " useDeltaR : ";
  switch(useDeltaR_) {
  case false : log << "false"; break;
  case true  : log << "true ";
  }
  log << "\n";
  log << " useMaxDist: ";
  switch(useMaxDist_) {
  case false : log << "false"; break;
  case true  : log << "true ";
  }
  log << "      maxDist: " << maxDist_ << "\n";
  log << " number of partons / jets: " << partons.size() << " / " << jets.size() << "\n";
  log << " number of available combinations: " << getNumberOfAvailableCombinations() << "\n";
  for(unsigned int comb = 0; comb < matching.size(); ++comb) {
    log << " -------------------------------------------- \n";
    log << " ind. of matched jets:";
    for(unsigned int part = 0; part < partons.size(); ++part)
      log << std::setw(4) << getMatchForParton(part, comb);
    log << "\n";
    log << " sumDeltaR             : " << getSumDeltaR(comb) << "\n";
    log << " sumDeltaPt / sumDeltaE: " << getSumDeltaPt(comb) << " / "  << getSumDeltaE(comb);
    log << "\n";
  }
  log << "++++++++++++++++++++++++++++++++++++++++++++++";
}