RecoB/src/MatchJet.cc

0001 #include <functional>
0002 #include <algorithm>
0003 #include <vector>
0004 #include <memory>
0005 #include <string>
0006
0007 #include <Math/GenVector/Cartesian3D.h>
0008 #include <Math/GenVector/VectorUtil.h>
0009
0010 #include "DQMOffline/RecoB/interface/MatchJet.h"
0011
0012 #include "Matching.h"
0013
0014 using namespace btag;
0015
0016 namespace {
0017   inline double sqr(double val) { return val * val; }
0018 }  // namespace
0019
0020 MatchJet::MatchJet(const edm::ParameterSet& pSet)
0021     : maxChi2(pSet.getParameter<double>("maxChi2")),
0022       sigmaDeltaR2(sqr(pSet.getParameter<double>("sigmaDeltaR"))),
0023       sigmaDeltaE2(sqr(pSet.getParameter<double>("sigmaDeltaE"))),
0024       threshold(1.0) {}
0025
0026 void MatchJet::matchCollections(const edm::RefToBaseVector<reco::Jet>& refJets_,
0027                                 const edm::RefToBaseVector<reco::Jet>& recJets_,
0028                                 const reco::JetCorrector* corrector) {
0029   refJetCorrector.setCorrector(corrector);
0030   recJetCorrector.setCorrector(corrector);
0031
0032   typedef ROOT::Math::Cartesian3D<double> Vector;
0033
0034   std::vector<Vector> corrRefJets;
0035   refJets.clear();
0036   for (edm::RefToBaseVector<reco::Jet>::const_iterator iter = refJets.begin(); iter != refJets_.end(); ++iter) {
0037     edm::RefToBase<reco::Jet> jetRef = *iter;
0038     reco::Jet jet = refJetCorrector(*jetRef);
0039     if (jet.energy() < threshold)
0040       continue;
0041
0042     corrRefJets.push_back(Vector(jet.px(), jet.py(), jet.pz()));
0043     refJets.push_back(jetRef);
0044   }
0045
0046   std::vector<Vector> corrRecJets;
0047   recJets.clear();
0048   for (edm::RefToBaseVector<reco::Jet>::const_iterator iter = recJets.begin(); iter != recJets_.end(); ++iter) {
0049     edm::RefToBase<reco::Jet> jetRec = *iter;
0050     reco::Jet jet = recJetCorrector(*jetRec);
0051     if (jet.energy() < threshold)
0052       continue;
0053
0054     corrRecJets.push_back(Vector(jet.px(), jet.py(), jet.pz()));
0055     recJets.push_back(jetRec);
0056   }
0057
0058   refToRec.clear();
0059   refToRec.resize(refJets.size(), -1);
0060
0061   recToRef.clear();
0062   recToRef.resize(recJets.size(), -1);
0063
0064   Matching<double> matching(corrRefJets, corrRecJets, [this](auto& v1, auto& v2) {
0065     double x = ROOT::Math::VectorUtil::DeltaR2(v1, v2) / this->sigmaDeltaR2 +
0066                sqr(2. * (v1.R() - v2.R()) / (v1.R() + v2.R())) / this->sigmaDeltaE2;
0067     // Possible debug output
0068     // std::cout << "xxx " << ROOT::Math::VectorUtil::DeltaPhi(v1, v2) << " " << (v1.Eta() - v2.Eta())
0069     // << " " << (v1.R() - v2.R()) / (v1.R() + v2.R()) << " " << x << std::endl;
0070     return x;
0071   });
0072
0073   typedef Matching<double>::Match Match;
0074
0075   const std::vector<Match>& matches = matching.match(std::less<double>(), [&](auto& c) { return c < this->maxChi2; });
0076   for (std::vector<Match>::const_iterator iter = matches.begin(); iter != matches.end(); ++iter) {
0077     refToRec[iter->index1] = iter->index2;
0078     recToRef[iter->index2] = iter->index1;
0079   }
0080 }
0081
0082 edm::RefToBase<reco::Jet> MatchJet::operator()(const edm::RefToBase<reco::Jet>& recJet) const {
0083   edm::RefToBase<reco::Jet> result;
0084   if (recJet.id() != recJets.id())
0085     return result;
0086
0087   for (unsigned int i = 0; i != recJets.size(); ++i) {
0088     if (recJets[i] == recJet) {
0089       int match = recToRef[i];
0090       if (match >= 0)
0091         result = refJets[match];
0092       break;
0093     }
0094   }
0095
0096   return result;
0097 }