Benchmark/interface/PFBenchmarkAlgo.h

0001 #ifndef RecoParticleFlow_Benchmark_PFBenchmarkAlgo_h
0002 #define RecoParticleFlow_Benchmark_PFBenchmarkAlgo_h
0003
0004 #include "FWCore/Utilities/interface/Exception.h"
0005 #include "DataFormats/Common/interface/OwnVector.h"
0006
0007 #include <cmath>
0008 #include <vector>
0009
0010 // Notes on template implementation:
0011 // - T, U are arbitrary types (must have et(), eta(), etc. defined)
0012 //   support for Candidate-derived objects is explicit
0013 // - Collection is a generic container class. Support for edm::OwnVector
0014 //   and std::vector is explicit.
0015
0016 class PFBenchmarkAlgo {
0017 public:
0018   // Calculate Delta-Et for the pair of Candidates (T - U)
0019   template <typename T, typename U>
0020   static double deltaEt(const T *, const U *);
0021
0022   // Calculate Delta-Eta for the pair of Candidates (T - U)
0023   template <typename T, typename U>
0024   static double deltaEta(const T *, const U *);
0025
0026   // Calculate Delta-Phi for the pair of Candidates (T - U)
0027   template <typename T, typename U>
0028   static double deltaPhi(const T *, const U *);
0029
0030   // Calculate Delta-R for the pair of Candidates
0031   template <typename T, typename U>
0032   static double deltaR(const T *, const U *);
0033
0034   // Match Candidate T to a Candidate in the Collection based on minimum Delta-R
0035   template <typename T, typename Collection>
0036   static const typename Collection::value_type *matchByDeltaR(const T *, const Collection *);
0037
0038   // Match Candidate T to a Candidate U in the Collection based on minimum Delta-Et
0039   template <typename T, typename Collection>
0040   static const typename Collection::value_type *matchByDeltaEt(const T *, const Collection *);
0041
0042   // Copy the input Collection (useful when sorting)
0043   template <typename T, typename Collection>
0044   static Collection copyCollection(const Collection *);
0045
0046   // Sort the U Candidates to the T Candidate based on minimum Delta-R
0047   template <typename T, typename Collection>
0048   static void sortByDeltaR(const T *, Collection &);
0049
0050   // Sort the U Candidates to the T Candidate based on minimum Delta-Et
0051   template <typename T, typename Collection>
0052   static void sortByDeltaEt(const T *, Collection &);
0053
0054   // Constrain the U Candidates to the T Candidate based on Delta-R to T
0055   template <typename T, typename Collection>
0056   static Collection findAllInCone(const T *, const Collection *, double);
0057
0058   // Constrain the U Candidates to the T Candidate based on Delta-Et to T
0059   template <typename T, typename Collection>
0060   static Collection findAllInEtWindow(const T *, const Collection *, double);
0061
0062 private:
0063   // std::vector sort helper function
0064   template <typename T, typename U, template <typename, typename> class Sorter>
0065   static void vector_sort(std::vector<T> &candidates, Sorter<T, U> S) {
0066     sort(candidates.begin(), candidates.end(), S);
0067   }
0068
0069   // std::vector push_back helper function
0070   template <typename T>
0071   static void vector_add(const T *c1, std::vector<T> &candidates) {
0072     candidates.push_back(*c1);
0073   }
0074
0075   // edm::OwnVector sort helper functions
0076   template <typename T, typename U, template <typename, typename> class Sorter>
0077   static void vector_sort(edm::OwnVector<T> &candidates, Sorter<T, U> S) {
0078     candidates.sort(S);
0079   }
0080
0081   // edm::OwnVector push_back helper function
0082   template <typename T>
0083   static void vector_add(const T *c1, edm::OwnVector<T> &candidates) {
0084     candidates.push_back((T *const)c1->clone());
0085   }
0086 };
0087
0088 // ========================================================================
0089 // implementation follows (required to be in header for templates)
0090 // ========================================================================
0091
0092 // Helper class for sorting U Collections by Delta-R to a Candidate T
0093 template <typename T, typename U>
0094 class deltaRSorter {
0095 public:
0096   deltaRSorter(const T *Ref) { cref = Ref; }
0097   bool operator()(const U &c1, const U &c2) const {
0098     return PFBenchmarkAlgo::deltaR(cref, &c1) < PFBenchmarkAlgo::deltaR(cref, &c2);
0099   }
0100
0101 private:
0102   const T *cref;
0103 };
0104
0105 // Helper class for sorting U Collections by Delta-Et to a Candidate T
0106 template <typename T, typename U>
0107 class deltaEtSorter {
0108 public:
0109   deltaEtSorter(const T *Ref) { cref = Ref; }
0110   bool operator()(const U &c1, const U &c2) const {
0111     return fabs(PFBenchmarkAlgo::deltaEt(cref, &c1)) < fabs(PFBenchmarkAlgo::deltaEt(cref, &c2));
0112   }
0113
0114 private:
0115   const T *cref;
0116 };
0117
0118 // Calculate Delta-Et for Candidates (T - U)
0119 template <typename T, typename U>
0120 double PFBenchmarkAlgo::deltaEt(const T *c1, const U *c2) {
0121   if (c1 == nullptr || c2 == nullptr)
0122     throw cms::Exception("Invalid Arg") << "attempted to calculate deltaEt for invalid Candidate(s)";
0123
0124   return c1->et() - c2->et();
0125 }
0126
0127 // Calculate Delta-Eta for Candidates (T - U)
0128 template <typename T, typename U>
0129 double PFBenchmarkAlgo::deltaEta(const T *c1, const U *c2) {
0130   if (c1 == nullptr || c2 == nullptr)
0131     throw cms::Exception("Invalid Arg") << "attempted to calculate deltaEta for invalid Candidate(s)";
0132
0133   return c1->eta() - c2->eta();
0134 }
0135
0136 // Calculate Delta-Phi for Candidates (T - U)
0137 template <typename T, typename U>
0138 double PFBenchmarkAlgo::deltaPhi(const T *c1, const U *c2) {
0139   if (c1 == nullptr || c2 == nullptr)
0140     throw cms::Exception("Invalid Arg") << "attempted to calculate deltaPhi for invalid Candidate(s)";
0141
0142   double phi1 = c1->phi();
0143   if (phi1 > M_PI)
0144     phi1 -= ceil((phi1 - M_PI) / (2 * M_PI)) * 2 * M_PI;
0145   if (phi1 <= -M_PI)
0146     phi1 += ceil((phi1 + M_PI) / (-2. * M_PI)) * 2. * M_PI;
0147
0148   double phi2 = c2->phi();
0149   if (phi2 > M_PI)
0150     phi2 -= ceil((phi2 - M_PI) / (2 * M_PI)) * 2 * M_PI;
0151   if (phi2 <= -M_PI)
0152     phi2 += ceil((phi2 + M_PI) / (-2. * M_PI)) * 2 * M_PI;
0153
0154   // alternative method:
0155   // while (phi > M_PI) phi -= 2 * M_PI;
0156   // while (phi <= - M_PI) phi += 2 * M_PI;
0157
0158   double deltaphi = -999.0;
0159   if (fabs(phi1 - phi2) < M_PI) {
0160     deltaphi = (phi1 - phi2);
0161   } else {
0162     if ((phi1 - phi2) > 0.0) {
0163       deltaphi = (2 * M_PI - fabs(phi1 - phi2));
0164     } else {
0165       deltaphi = -(2 * M_PI - fabs(phi1 - phi2));
0166     }
0167   }
0168   return deltaphi;
0169   //return ( (fabs(phi1 - phi2)<M_PI)?(phi1-phi2):(2*M_PI - fabs(phi1 - phi2) ) ); // FL: wrong
0170 }
0171
0172 // Calculate Delta-R for Candidates
0173 template <typename T, typename U>
0174 double PFBenchmarkAlgo::deltaR(const T *c1, const U *c2) {
0175   if (c1 == nullptr || c2 == nullptr)
0176     throw cms::Exception("Invalid Arg") << "attempted to calculate deltaR for invalid Candidate(s)";
0177
0178   return sqrt(std::pow(deltaPhi(c1, c2), 2) + std::pow(deltaEta(c1, c2), 2));
0179 }
0180
0181 // Match Candidate T to a Candidate in the Collection based on minimum Delta-R
0182 template <typename T, typename Collection>
0183 const typename Collection::value_type *PFBenchmarkAlgo::matchByDeltaR(const T *c1, const Collection *candidates) {
0184   typedef typename Collection::value_type U;
0185
0186   // Try to verify the validity of the Candidate and Collection
0187   if (!c1)
0188     throw cms::Exception("Invalid Arg") << "attempted to match invalid Candidate";
0189   if (!candidates)
0190     throw cms::Exception("Invalid Arg") << "attempted to match to invalid Collection";
0191
0192   double minDeltaR = 9999.;
0193   const U *match = nullptr;
0194
0195   // Loop Over the Candidates...
0196   for (unsigned int i = 0; i < candidates->size(); i++) {
0197     const U *c2 = &(*candidates)[i];
0198     if (!c2)
0199       throw cms::Exception("Invalid Arg") << "attempted to match to invalid Candidate";
0200
0201     // Find Minimal Delta-R
0202     double dR = deltaR(c1, c2);
0203     if (dR <= minDeltaR) {
0204       match = c2;
0205       minDeltaR = dR;
0206     }
0207   }
0208
0209   // Return the Candidate with the smallest dR
0210   return match;
0211 }
0212
0213 // Match Candidate T to a Candidate U in the Collection based on minimum Delta-Et
0214 template <typename T, typename Collection>
0215 const typename Collection::value_type *PFBenchmarkAlgo::matchByDeltaEt(const T *c1, const Collection *candidates) {
0216   typedef typename Collection::value_type U;
0217
0218   // Try to verify the validity of the Candidate and Collection
0219   if (!c1)
0220     throw cms::Exception("Invalid Arg") << "attempted to match invalid Candidate";
0221   if (!candidates)
0222     throw cms::Exception("Invalid Arg") << "attempted to match to invalid Collection";
0223
0224   double minDeltaEt = 9999.;
0225   const U *match = NULL;
0226
0227   // Loop Over the Candidates...
0228   for (unsigned int i = 0; i < candidates->size(); i++) {
0229     const T *c2 = &(*candidates)[i];
0230     if (!c2)
0231       throw cms::Exception("Invalid Arg") << "attempted to match to invalid Candidate";
0232
0233     // Find Minimal Delta-R
0234     double dEt = fabs(deltaEt(c1, c2));
0235     if (dEt <= minDeltaEt) {
0236       match = c2;
0237       minDeltaEt = dEt;
0238     }
0239   }
0240
0241   // Return the Candidate with the smallest dR
0242   return match;
0243 }
0244
0245 // Copy the Collection (useful when sorting)
0246 template <typename T, typename Collection>
0247 Collection PFBenchmarkAlgo::copyCollection(const Collection *candidates) {
0248   typedef typename Collection::value_type U;
0249
0250   // Try to verify the validity of the Collection
0251   if (!candidates)
0252     throw cms::Exception("Invalid Arg") << "attempted to copy invalid Collection";
0253
0254   Collection copy;
0255
0256   for (unsigned int i = 0; i < candidates->size(); i++)
0257     vector_add(&(*candidates)[i], copy);
0258
0259   return copy;
0260 }
0261
0262 // Sort the U Candidates to the Candidate T based on minimum Delta-R
0263 template <typename T, typename Collection>
0264 void PFBenchmarkAlgo::sortByDeltaR(const T *c1, Collection &candidates) {
0265   typedef typename Collection::value_type U;
0266
0267   // Try to verify the validity of Candidate and Collection
0268   if (!c1)
0269     throw cms::Exception("Invalid Arg") << "attempted to sort by invalid Candidate";
0270   if (!candidates)
0271     throw cms::Exception("Invalid Arg") << "attempted to sort invalid Candidates";
0272
0273   // Sort the collection
0274   vector_sort(candidates, deltaRSorter<T, U>(c1));
0275 }
0276
0277 // Sort the U Candidates to the Candidate T based on minimum Delta-Et
0278 template <typename T, typename Collection>
0279 void PFBenchmarkAlgo::sortByDeltaEt(const T *c1, Collection &candidates) {
0280   typedef typename Collection::value_type U;
0281
0282   // Try to verify the validity of Candidate and Collection
0283   if (!c1)
0284     throw cms::Exception("Invalid Arg") << "attempted to sort by invalid Candidate";
0285   if (!candidates)
0286     throw cms::Exception("Invalid Arg") << "attempted to sort invalid Candidates";
0287
0288   // Sort the collection
0289   vector_sort(candidates, deltaEtSorter<T, U>(c1));
0290 }
0291
0292 // Constrain the U Candidates to the T Candidate based on Delta-R to T
0293 template <typename T, typename Collection>
0294 Collection PFBenchmarkAlgo::findAllInCone(const T *c1, const Collection *candidates, double ConeSize) {
0295   typedef typename Collection::value_type U;
0296
0297   // Try to verify the validity of Candidate and the Collection
0298   if (!c1)
0299     throw cms::Exception("Invalid Arg") << "attempted to constrain to invalid Candidate";
0300   if (!candidates)
0301     throw cms::Exception("Invalid Arg") << "attempted to constrain invalid Collection";
0302   if (ConeSize <= 0)
0303     throw cms::Exception("Invalid Arg") << "zero or negative cone size specified";
0304
0305   Collection constrained;
0306
0307   for (unsigned int i = 0; i < candidates->size(); i++) {
0308     const U *c2 = &(*candidates)[i];
0309
0310     // Add in-cone Candidates to the new Collection
0311     double dR = deltaR(c1, c2);
0312     if (dR < ConeSize)
0313       vector_add(c2, constrained);
0314   }
0315
0316   // Sort and return Collection
0317   sortByDeltaR(c1, constrained);
0318   return constrained;
0319 }
0320
0321 // Constrain the U Candidates to the T Candidate based on Delta-Et to T
0322 template <typename T, typename Collection>
0323 Collection PFBenchmarkAlgo::findAllInEtWindow(const T *c1, const Collection *candidates, double EtWindow) {
0324   typedef typename Collection::value_type U;
0325
0326   // Try to verify the validity of Candidate and the Collection
0327   if (!c1)
0328     throw cms::Exception("Invalid Arg") << "attempted to constrain to invalid Candidate";
0329   if (!candidates)
0330     throw cms::Exception("Invalid Arg") << "attempted to constrain invalid Collection";
0331   if (EtWindow <= 0)
0332     throw cms::Exception("Invalid Arg") << "zero or negative cone size specified";
0333
0334   Collection constrained;
0335
0336   //CandidateCollection::const_iterator candidate;
0337   //for (candidate = candidates->begin(); candidate != candidates->end(); candidate++) {
0338   for (unsigned int i = 0; i < candidates->size(); i++) {
0339     const U *c2 = &(*candidates)[i];
0340
0341     // Add in-cone Candidates to the new Collection
0342     double dEt = fabs(deltaEt(c1, c2));
0343     if (dEt < EtWindow)
0344       vector_add(c2, constrained);
0345   }
0346
0347   // Sort and return Collection
0348   sortByDeltaEt(c1, constrained);
0349   return constrained;
0350 }
0351
0352 #endif  // RecoParticleFlow_Benchmark_PFBenchmarkAlgo_h