PFBenchmarkAlgo.h

Go to the documentation of this file.

#ifndef RecoParticleFlow_Benchmark_PFBenchmarkAlgo_h
#define RecoParticleFlow_Benchmark_PFBenchmarkAlgo_h

#include "FWCore/Utilities/interface/Exception.h"
#include "DataFormats/Common/interface/OwnVector.h"

#include <cmath>
#include <vector>

// Notes on template implementation:
// - T, U are arbitrary types (must have et(), eta(), etc. defined)
//   support for Candidate-derived objects is explicit
// - Collection is a generic container class. Support for edm::OwnVector
//   and std::vector is explicit.

class PFBenchmarkAlgo {
public:
  // Calculate Delta-Et for the pair of Candidates (T - U)
  template <typename T, typename U>
  static double deltaEt(const T *, const U *);

  // Calculate Delta-Eta for the pair of Candidates (T - U)
  template <typename T, typename U>
  static double deltaEta(const T *, const U *);

  // Calculate Delta-Phi for the pair of Candidates (T - U)
  template <typename T, typename U>
  static double deltaPhi(const T *, const U *);

  // Calculate Delta-R for the pair of Candidates
  template <typename T, typename U>
  static double deltaR(const T *, const U *);

  // Match Candidate T to a Candidate in the Collection based on minimum Delta-R
  template <typename T, typename Collection>
  static const typename Collection::value_type *matchByDeltaR(const T *, const Collection *);

  // Match Candidate T to a Candidate U in the Collection based on minimum Delta-Et
  template <typename T, typename Collection>
  static const typename Collection::value_type *matchByDeltaEt(const T *, const Collection *);

  // Copy the input Collection (useful when sorting)
  template <typename T, typename Collection>
  static Collection copyCollection(const Collection *);

  // Sort the U Candidates to the T Candidate based on minimum Delta-R
  template <typename T, typename Collection>
  static void sortByDeltaR(const T *, Collection &);

  // Sort the U Candidates to the T Candidate based on minimum Delta-Et
  template <typename T, typename Collection>
  static void sortByDeltaEt(const T *, Collection &);

  // Constrain the U Candidates to the T Candidate based on Delta-R to T
  template <typename T, typename Collection>
  static Collection findAllInCone(const T *, const Collection *, double);

  // Constrain the U Candidates to the T Candidate based on Delta-Et to T
  template <typename T, typename Collection>
  static Collection findAllInEtWindow(const T *, const Collection *, double);

private:
  // std::vector sort helper function
  template <typename T, typename U, template <typename, typename> class Sorter>
  static void vector_sort(std::vector<T> &candidates, Sorter<T, U> S) {
    sort(candidates.begin(), candidates.end(), S);
  }

  // std::vector push_back helper function
  template <typename T>
  static void vector_add(const T *c1, std::vector<T> &candidates) {
    candidates.push_back(*c1);
  }

  // edm::OwnVector sort helper functions
  template <typename T, typename U, template <typename, typename> class Sorter>
  static void vector_sort(edm::OwnVector<T> &candidates, Sorter<T, U> S) {
    candidates.sort(S);
  }

  // edm::OwnVector push_back helper function
  template <typename T>
  static void vector_add(const T *c1, edm::OwnVector<T> &candidates) {
    candidates.push_back((T *const)c1->clone());
  }
};

// ========================================================================
// implementation follows (required to be in header for templates)
// ========================================================================

// Helper class for sorting U Collections by Delta-R to a Candidate T
template <typename T, typename U>
class deltaRSorter {
public:
  deltaRSorter(const T *Ref) { cref = Ref; }
  bool operator()(const U &c1, const U &c2) const {
    return PFBenchmarkAlgo::deltaR(cref, &c1) < PFBenchmarkAlgo::deltaR(cref, &c2);
  }

private:
  const T *cref;
};

// Helper class for sorting U Collections by Delta-Et to a Candidate T
template <typename T, typename U>
class deltaEtSorter {
public:
  deltaEtSorter(const T *Ref) { cref = Ref; }
  bool operator()(const U &c1, const U &c2) const {
    return fabs(PFBenchmarkAlgo::deltaEt(cref, &c1)) < fabs(PFBenchmarkAlgo::deltaEt(cref, &c2));
  }

private:
  const T *cref;
};

// Calculate Delta-Et for Candidates (T - U)
template <typename T, typename U>
double PFBenchmarkAlgo::deltaEt(const T *c1, const U *c2) {
  if (c1 == nullptr || c2 == nullptr)
    throw cms::Exception("Invalid Arg") << "attempted to calculate deltaEt for invalid Candidate(s)";

  return c1->et() - c2->et();
}

// Calculate Delta-Eta for Candidates (T - U)
template <typename T, typename U>
double PFBenchmarkAlgo::deltaEta(const T *c1, const U *c2) {
  if (c1 == nullptr || c2 == nullptr)
    throw cms::Exception("Invalid Arg") << "attempted to calculate deltaEta for invalid Candidate(s)";

  return c1->eta() - c2->eta();
}

// Calculate Delta-Phi for Candidates (T - U)
template <typename T, typename U>
double PFBenchmarkAlgo::deltaPhi(const T *c1, const U *c2) {
  if (c1 == nullptr || c2 == nullptr)
    throw cms::Exception("Invalid Arg") << "attempted to calculate deltaPhi for invalid Candidate(s)";

  double phi1 = c1->phi();
  if (phi1 > M_PI)
    phi1 -= ceil((phi1 - M_PI) / (2 * M_PI)) * 2 * M_PI;
  if (phi1 <= -M_PI)
    phi1 += ceil((phi1 + M_PI) / (-2. * M_PI)) * 2. * M_PI;

  double phi2 = c2->phi();
  if (phi2 > M_PI)
    phi2 -= ceil((phi2 - M_PI) / (2 * M_PI)) * 2 * M_PI;
  if (phi2 <= -M_PI)
    phi2 += ceil((phi2 + M_PI) / (-2. * M_PI)) * 2 * M_PI;

  // alternative method:
  // while (phi > M_PI) phi -= 2 * M_PI;
  // while (phi <= - M_PI) phi += 2 * M_PI;

  double deltaphi = -999.0;
  if (fabs(phi1 - phi2) < M_PI) {
    deltaphi = (phi1 - phi2);
  } else {
    if ((phi1 - phi2) > 0.0) {
      deltaphi = (2 * M_PI - fabs(phi1 - phi2));
    } else {
      deltaphi = -(2 * M_PI - fabs(phi1 - phi2));
    }
  }
  return deltaphi;
  //return ( (fabs(phi1 - phi2)<M_PI)?(phi1-phi2):(2*M_PI - fabs(phi1 - phi2) ) ); // FL: wrong
}

// Calculate Delta-R for Candidates
template <typename T, typename U>
double PFBenchmarkAlgo::deltaR(const T *c1, const U *c2) {
  if (c1 == nullptr || c2 == nullptr)
    throw cms::Exception("Invalid Arg") << "attempted to calculate deltaR for invalid Candidate(s)";

  return sqrt(std::pow(deltaPhi(c1, c2), 2) + std::pow(deltaEta(c1, c2), 2));
}

// Match Candidate T to a Candidate in the Collection based on minimum Delta-R
template <typename T, typename Collection>
const typename Collection::value_type *PFBenchmarkAlgo::matchByDeltaR(const T *c1, const Collection *candidates) {
  typedef typename Collection::value_type U;

  // Try to verify the validity of the Candidate and Collection
  if (!c1)
    throw cms::Exception("Invalid Arg") << "attempted to match invalid Candidate";
  if (!candidates)
    throw cms::Exception("Invalid Arg") << "attempted to match to invalid Collection";

  double minDeltaR = 9999.;
  const U *match = nullptr;

  // Loop Over the Candidates...
  for (unsigned int i = 0; i < candidates->size(); i++) {
    const U *c2 = &(*candidates)[i];
    if (!c2)
      throw cms::Exception("Invalid Arg") << "attempted to match to invalid Candidate";

    // Find Minimal Delta-R
    double dR = deltaR(c1, c2);
    if (dR <= minDeltaR) {
      match = c2;
      minDeltaR = dR;
    }
  }

  // Return the Candidate with the smallest dR
  return match;
}

// Match Candidate T to a Candidate U in the Collection based on minimum Delta-Et
template <typename T, typename Collection>
const typename Collection::value_type *PFBenchmarkAlgo::matchByDeltaEt(const T *c1, const Collection *candidates) {
  typedef typename Collection::value_type U;

  // Try to verify the validity of the Candidate and Collection
  if (!c1)
    throw cms::Exception("Invalid Arg") << "attempted to match invalid Candidate";
  if (!candidates)
    throw cms::Exception("Invalid Arg") << "attempted to match to invalid Collection";

  double minDeltaEt = 9999.;
  const U *match = NULL;

  // Loop Over the Candidates...
  for (unsigned int i = 0; i < candidates->size(); i++) {
    const T *c2 = &(*candidates)[i];
    if (!c2)
      throw cms::Exception("Invalid Arg") << "attempted to match to invalid Candidate";

    // Find Minimal Delta-R
    double dEt = fabs(deltaEt(c1, c2));
    if (dEt <= minDeltaEt) {
      match = c2;
      minDeltaEt = dEt;
    }
  }

  // Return the Candidate with the smallest dR
  return match;
}

// Copy the Collection (useful when sorting)
template <typename T, typename Collection>
Collection PFBenchmarkAlgo::copyCollection(const Collection *candidates) {
  typedef typename Collection::value_type U;

  // Try to verify the validity of the Collection
  if (!candidates)
    throw cms::Exception("Invalid Arg") << "attempted to copy invalid Collection";

  Collection copy;

  for (unsigned int i = 0; i < candidates->size(); i++)
    vector_add(&(*candidates)[i], copy);

  return copy;
}

// Sort the U Candidates to the Candidate T based on minimum Delta-R
template <typename T, typename Collection>
void PFBenchmarkAlgo::sortByDeltaR(const T *c1, Collection &candidates) {
  typedef typename Collection::value_type U;

  // Try to verify the validity of Candidate and Collection
  if (!c1)
    throw cms::Exception("Invalid Arg") << "attempted to sort by invalid Candidate";
  if (!candidates)
    throw cms::Exception("Invalid Arg") << "attempted to sort invalid Candidates";

  // Sort the collection
  vector_sort(candidates, deltaRSorter<T, U>(c1));
}

// Sort the U Candidates to the Candidate T based on minimum Delta-Et
template <typename T, typename Collection>
void PFBenchmarkAlgo::sortByDeltaEt(const T *c1, Collection &candidates) {
  typedef typename Collection::value_type U;

  // Try to verify the validity of Candidate and Collection
  if (!c1)
    throw cms::Exception("Invalid Arg") << "attempted to sort by invalid Candidate";
  if (!candidates)
    throw cms::Exception("Invalid Arg") << "attempted to sort invalid Candidates";

  // Sort the collection
  vector_sort(candidates, deltaEtSorter<T, U>(c1));
}

// Constrain the U Candidates to the T Candidate based on Delta-R to T
template <typename T, typename Collection>
Collection PFBenchmarkAlgo::findAllInCone(const T *c1, const Collection *candidates, double ConeSize) {
  typedef typename Collection::value_type U;

  // Try to verify the validity of Candidate and the Collection
  if (!c1)
    throw cms::Exception("Invalid Arg") << "attempted to constrain to invalid Candidate";
  if (!candidates)
    throw cms::Exception("Invalid Arg") << "attempted to constrain invalid Collection";
  if (ConeSize <= 0)
    throw cms::Exception("Invalid Arg") << "zero or negative cone size specified";

  Collection constrained;

  for (unsigned int i = 0; i < candidates->size(); i++) {
    const U *c2 = &(*candidates)[i];

    // Add in-cone Candidates to the new Collection
    double dR = deltaR(c1, c2);
    if (dR < ConeSize)
      vector_add(c2, constrained);
  }

  // Sort and return Collection
  sortByDeltaR(c1, constrained);
  return constrained;
}

// Constrain the U Candidates to the T Candidate based on Delta-Et to T
template <typename T, typename Collection>
Collection PFBenchmarkAlgo::findAllInEtWindow(const T *c1, const Collection *candidates, double EtWindow) {
  typedef typename Collection::value_type U;

  // Try to verify the validity of Candidate and the Collection
  if (!c1)
    throw cms::Exception("Invalid Arg") << "attempted to constrain to invalid Candidate";
  if (!candidates)
    throw cms::Exception("Invalid Arg") << "attempted to constrain invalid Collection";
  if (EtWindow <= 0)
    throw cms::Exception("Invalid Arg") << "zero or negative cone size specified";

  Collection constrained;

  //CandidateCollection::const_iterator candidate;
  //for (candidate = candidates->begin(); candidate != candidates->end(); candidate++) {
  for (unsigned int i = 0; i < candidates->size(); i++) {
    const U *c2 = &(*candidates)[i];

    // Add in-cone Candidates to the new Collection
    double dEt = fabs(deltaEt(c1, c2));
    if (dEt < EtWindow)
      vector_add(c2, constrained);
  }

  // Sort and return Collection
  sortByDeltaEt(c1, constrained);
  return constrained;
}

#endif  // RecoParticleFlow_Benchmark_PFBenchmarkAlgo_h