UtilAlgos/interface/PhysObjectMatcher.h

0001 #ifndef UtilAlgos_PhysObjectMatcher_h
0002 #define UtilAlgos_PhysObjectMatcher_h
0003 /* \class PhysObjectMatcher.
0004  *
0005  * Extended version of reco::CandMatcher.
0006  * Tries to match elements from collection 1 to collection 2 with optional
0007  * resolution of ambiguities. Uses three helper classes for
0008  *  (1) the basic selection of the match (e.g. pdgId, charge, ..);
0009  *  (2) a distance measure (e.g. deltaR);
0010  *  (3) the ranking of several matches.
0011  *
0012  */
0013 #include "FWCore/Framework/interface/stream/EDProducer.h"
0014 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0015 #include "FWCore/Utilities/interface/InputTag.h"
0016 #include "CommonTools/UtilAlgos/interface/DeltaR.h"
0017 #include "CommonTools/UtilAlgos/interface/MasterCollectionHelper.h"
0018 #include "FWCore/Framework/interface/Event.h"
0019 #include "DataFormats/Common/interface/Handle.h"
0020 #include "DataFormats/Common/interface/Association.h"
0021 #include "CommonTools/UtilAlgos/interface/MatchByDR.h"
0022
0023 // #include <iostream>
0024
0025 namespace reco {
0026
0027   namespace helper {
0028     /// Default class for ranking matches: sorting by smaller distance.
0029     template <typename D, typename C1, typename C2>
0030     class LessByMatchDistance {
0031     public:
0032       LessByMatchDistance(const edm::ParameterSet& cfg, const C1& c1, const C2& c2)
0033           : distance_(reco::modules::make<D>(cfg)), c1_(c1), c2_(c2) {}
0034       bool operator()(const std::pair<size_t, size_t>& p1, const std::pair<size_t, size_t>& p2) const {
0035         return distance_(c1_[p1.first], c2_[p1.second]) < distance_(c1_[p2.first], c2_[p2.second]);
0036       }
0037
0038     private:
0039       D distance_;
0040       const C1& c1_;
0041       const C2& c2_;
0042     };
0043   }  // namespace helper
0044
0045   // Template arguments:
0046   // C1 .. candidate collection to be matched
0047   // C2 .. target of the match (typically MC)
0048   // S ... match (pre-)selector
0049   // D ... match (typically cut on some distance)
0050   //         default: deltaR cut
0051   // Q ... ranking of matches
0052   //         default: by smaller deltaR
0053   template <typename C1,
0054             typename C2,
0055             typename S,
0056             typename D = reco::MatchByDR<typename C1::value_type, typename C2::value_type>,
0057             typename Q = helper::LessByMatchDistance<DeltaR<typename C1::value_type, typename C2::value_type>, C1, C2> >
0058   class PhysObjectMatcher : public edm::stream::EDProducer<> {
0059   public:
0060     PhysObjectMatcher(const edm::ParameterSet& cfg);
0061     ~PhysObjectMatcher() override;
0062
0063   private:
0064     typedef typename C1::value_type T1;
0065     typedef typename C2::value_type T2;
0066     typedef edm::Association<C2> MatchMap;
0067     typedef std::pair<size_t, size_t> IndexPair;
0068     typedef std::vector<IndexPair> MatchContainer;
0069     void produce(edm::Event&, const edm::EventSetup&) override;
0070     edm::ParameterSet config_;
0071     edm::EDGetTokenT<C1> srcToken_;
0072     edm::EDGetTokenT<C2> matchedToken_;
0073     bool resolveAmbiguities_;     // resolve ambiguities after
0074                                   //   first pass?
0075     bool resolveByMatchQuality_;  // resolve by (global) quality
0076                                   //   of match (otherwise: by order
0077                                   //   of test candidates)
0078     bool select(const T1& c1, const T2& c2) const { return select_(c1, c2); }
0079     S select_;
0080     D distance_;
0081     //     DeltaR<typename C1::value_type, typename C2::value_type> testDR_;
0082   };
0083
0084   template <typename C1, typename C2, typename S, typename D, typename Q>
0085   PhysObjectMatcher<C1, C2, S, D, Q>::PhysObjectMatcher(const edm::ParameterSet& cfg)
0086       : config_(cfg),
0087         srcToken_(consumes<C1>(cfg.template getParameter<edm::InputTag>("src"))),
0088         matchedToken_(consumes<C2>(cfg.template getParameter<edm::InputTag>("matched"))),
0089         resolveAmbiguities_(cfg.template getParameter<bool>("resolveAmbiguities")),
0090         resolveByMatchQuality_(cfg.template getParameter<bool>("resolveByMatchQuality")),
0091         select_(reco::modules::make<S>(cfg)),
0092         distance_(reco::modules::make<D>(cfg)) {
0093     // definition of the product
0094     produces<MatchMap>();
0095     // set resolveByMatchQuality only if ambiguities are to be resolved
0096     resolveByMatchQuality_ = resolveByMatchQuality_ && resolveAmbiguities_;
0097   }
0098
0099   template <typename C1, typename C2, typename S, typename D, typename Q>
0100   PhysObjectMatcher<C1, C2, S, D, Q>::~PhysObjectMatcher() {}
0101
0102   template <typename C1, typename C2, typename S, typename D, typename Q>
0103   void PhysObjectMatcher<C1, C2, S, D, Q>::produce(edm::Event& evt, const edm::EventSetup&) {
0104     using namespace edm;
0105     using namespace std;
0106     typedef std::pair<size_t, size_t> IndexPair;
0107     typedef std::vector<IndexPair> MatchContainer;
0108     // get collections from event
0109     Handle<C2> matched;
0110     evt.getByToken(matchedToken_, matched);
0111     Handle<C1> cands;
0112     evt.getByToken(srcToken_, cands);
0113     // create product
0114     unique_ptr<MatchMap> matchMap(new MatchMap(matched));
0115     size_t size = cands->size();
0116     if (size != 0) {
0117       //
0118       // create helpers
0119       //
0120       Q comparator(config_, *cands, *matched);
0121       typename MatchMap::Filler filler(*matchMap);
0122       ::helper::MasterCollection<C1> master(cands, evt);
0123       vector<int> indices(master.size(), -1);      // result: indices in target collection
0124       vector<bool> mLock(matched->size(), false);  // locks in target collection
0125       MatchContainer matchPairs;                   // container of matched pairs
0126       // loop over candidates
0127       for (size_t c = 0; c != size; ++c) {
0128         const T1& cand = (*cands)[c];
0129         // no global comparison of match quality -> reset the container for each candidate
0130         if (!resolveByMatchQuality_)
0131           matchPairs.clear();
0132         // loop over target collection
0133         for (size_t m = 0; m != matched->size(); ++m) {
0134           const T2& match = (*matched)[m];
0135           // check lock and preselection
0136           if (!mLock[m] && select(cand, match)) {
0137             //          double dist = testDR_(cand,match);
0138             //          cout << "dist between c = " << c << " and m = "
0139             //           << m << " is " << dist << " at pts of "
0140             //           << cand.pt() << " " << match.pt() << endl;
0141             // matching requirement fulfilled -> store pair of indices
0142             if (distance_(cand, match))
0143               matchPairs.push_back(make_pair(c, m));
0144           }
0145         }
0146         // if match(es) found and no global ambiguity resolution requested
0147         if (!matchPairs.empty() && !resolveByMatchQuality_) {
0148           // look for and store best match
0149           size_t idx = master.index(c);
0150           assert(idx < indices.size());
0151           size_t index = min_element(matchPairs.begin(), matchPairs.end(), comparator)->second;
0152           indices[idx] = index;
0153           // if ambiguity resolution by order of (reco) candidates:
0154           //   lock element in target collection
0155           if (resolveAmbiguities_)
0156             mLock[index] = true;
0157           //      {
0158           //        MatchContainer::const_iterator i = min_element(matchPairs.begin(), matchPairs.end(), comparator);
0159           //        cout << "smallest distance for c = " << c << " is "
0160           //         << testDR_((*cands)[(*i).first],
0161           //                (*matched)[(*i).second]) << endl;
0162           //      }
0163         }
0164       }
0165       // ambiguity resolution by global match quality (if requested)
0166       if (resolveByMatchQuality_) {
0167         // sort container of all matches by quality
0168         sort(matchPairs.begin(), matchPairs.end(), comparator);
0169         vector<bool> cLock(master.size(), false);
0170         // loop over sorted container
0171         for (MatchContainer::const_iterator i = matchPairs.begin(); i != matchPairs.end(); ++i) {
0172           size_t c = (*i).first;
0173           size_t m = (*i).second;
0174           //      cout << "rel dp = " << ((*cands)[c].pt()-(*matched)[m].pt())/(*matched)[m].pt() << endl;
0175           // accept only pairs without any lock
0176           if (mLock[m] || cLock[c])
0177             continue;
0178           // store index to target collection and lock the two items
0179           size_t idx = master.index(c);
0180           assert(idx < indices.size());
0181           indices[idx] = m;
0182           mLock[m] = true;
0183           cLock[c] = true;
0184         }
0185       }
0186       filler.insert(master.get(), indices.begin(), indices.end());
0187       filler.fill();
0188     }
0189     evt.put(std::move(matchMap));
0190   }
0191
0192 }  // namespace reco
0193
0194 #endif