File indexing completed on 2024-04-06 12:31:25
0001 #include "TopQuarkAnalysis/TopTools/interface/JetPartonMatching.h"
0002 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0003 #include <Math/VectorUtil.h>
0004
0005 JetPartonMatching::JetPartonMatching(const std::vector<const reco::Candidate*>& p,
0006 const std::vector<reco::GenJet>& j,
0007 const int algorithm = totalMinDist,
0008 const bool useMaxDist = true,
0009 const bool useDeltaR = true,
0010 const double maxDist = 0.3)
0011 : partons(p), algorithm_(algorithm), useMaxDist_(useMaxDist), useDeltaR_(useDeltaR), maxDist_(maxDist) {
0012 std::vector<const reco::Candidate*> js;
0013 for (unsigned int i = 0; i < j.size(); ++i)
0014 js.push_back(&(j[i]));
0015 jets = js;
0016 calculate();
0017 }
0018
0019 JetPartonMatching::JetPartonMatching(const std::vector<const reco::Candidate*>& p,
0020 const std::vector<reco::CaloJet>& j,
0021 const int algorithm = totalMinDist,
0022 const bool useMaxDist = true,
0023 const bool useDeltaR = true,
0024 const double maxDist = 0.3)
0025 : partons(p), algorithm_(algorithm), useMaxDist_(useMaxDist), useDeltaR_(useDeltaR), maxDist_(maxDist) {
0026 std::vector<const reco::Candidate*> js;
0027 for (unsigned int i = 0; i < j.size(); ++i)
0028 js.push_back(&(j[i]));
0029 jets = js;
0030 calculate();
0031 }
0032
0033 JetPartonMatching::JetPartonMatching(const std::vector<const reco::Candidate*>& p,
0034 const std::vector<pat::Jet>& j,
0035 const int algorithm = totalMinDist,
0036 const bool useMaxDist = true,
0037 const bool useDeltaR = true,
0038 const double maxDist = 0.3)
0039 : partons(p), algorithm_(algorithm), useMaxDist_(useMaxDist), useDeltaR_(useDeltaR), maxDist_(maxDist) {
0040 std::vector<const reco::Candidate*> js;
0041 for (unsigned int i = 0; i < j.size(); ++i)
0042 js.push_back(&(j[i]));
0043 jets = js;
0044 calculate();
0045 }
0046
0047 JetPartonMatching::JetPartonMatching(const std::vector<const reco::Candidate*>& p,
0048 const std::vector<const reco::Candidate*>& j,
0049 const int algorithm = totalMinDist,
0050 const bool useMaxDist = true,
0051 const bool useDeltaR = true,
0052 const double maxDist = 0.3)
0053 : partons(p), jets(j), algorithm_(algorithm), useMaxDist_(useMaxDist), useDeltaR_(useDeltaR), maxDist_(maxDist) {
0054 calculate();
0055 }
0056
0057 void JetPartonMatching::calculate() {
0058
0059
0060 if (algorithm_ == unambiguousOnly)
0061 useMaxDist_ = true;
0062
0063
0064
0065
0066
0067 bool emptyParton = false;
0068 for (unsigned int ip = 0; ip < partons.size(); ++ip) {
0069 if (partons[ip]->pdgId() == 0) {
0070 emptyParton = true;
0071 break;
0072 }
0073 }
0074
0075
0076
0077
0078 if (jets.empty() || emptyParton) {
0079 MatchingCollection dummyMatch;
0080 for (unsigned int ip = 0; ip < partons.size(); ++ip)
0081 dummyMatch.push_back(std::make_pair(ip, -1));
0082 matching.push_back(dummyMatch);
0083 } else {
0084 switch (algorithm_) {
0085 case totalMinDist:
0086 matchingTotalMinDist();
0087 break;
0088
0089 case minSumDist:
0090 matchingMinSumDist();
0091 break;
0092
0093 case ptOrderedMinDist:
0094 matchingPtOrderedMinDist();
0095 break;
0096
0097 case unambiguousOnly:
0098 matchingUnambiguousOnly();
0099 break;
0100
0101 default:
0102 matchingMinSumDist();
0103 }
0104 }
0105
0106 numberOfUnmatchedPartons.clear();
0107 sumDeltaE.clear();
0108 sumDeltaPt.clear();
0109 sumDeltaR.clear();
0110 for (unsigned int comb = 0; comb < matching.size(); ++comb) {
0111 MatchingCollection match = matching[comb];
0112 std::sort(match.begin(), match.end());
0113 matching[comb] = match;
0114
0115 int nUnmatchedPartons = partons.size();
0116 for (unsigned int part = 0; part < partons.size(); ++part)
0117 if (getMatchForParton(part, comb) >= 0)
0118 --nUnmatchedPartons;
0119
0120 double sumDE = -999.;
0121 double sumDPt = -999.;
0122 double sumDR = -999.;
0123 if (nUnmatchedPartons == 0) {
0124 sumDE = 0;
0125 sumDPt = 0;
0126 sumDR = 0;
0127 for (unsigned int i = 0; i < match.size(); ++i) {
0128 sumDE += fabs(partons[match[i].first]->energy() - jets[match[i].second]->energy());
0129 sumDPt += fabs(partons[match[i].first]->pt() - jets[match[i].second]->pt());
0130 sumDR += distance(partons[match[i].first]->p4(), jets[match[i].second]->p4());
0131 }
0132 }
0133
0134 numberOfUnmatchedPartons.push_back(nUnmatchedPartons);
0135 sumDeltaE.push_back(sumDE);
0136 sumDeltaPt.push_back(sumDPt);
0137 sumDeltaR.push_back(sumDR);
0138 }
0139 }
0140
0141 double JetPartonMatching::distance(const math::XYZTLorentzVector& v1, const math::XYZTLorentzVector& v2) {
0142
0143
0144 if (useDeltaR_)
0145 return ROOT::Math::VectorUtil::DeltaR(v1, v2);
0146 return ROOT::Math::VectorUtil::Angle(v1, v2);
0147 }
0148
0149 void JetPartonMatching::matchingTotalMinDist() {
0150
0151
0152
0153
0154
0155
0156 std::vector<std::pair<double, unsigned int> > distances;
0157 for (unsigned int ip = 0; ip < partons.size(); ++ip) {
0158 for (unsigned int ij = 0; ij < jets.size(); ++ij) {
0159 double dist = distance(jets[ij]->p4(), partons[ip]->p4());
0160 distances.push_back(std::pair<double, unsigned int>(dist, ip * jets.size() + ij));
0161 }
0162 }
0163 std::sort(distances.begin(), distances.end());
0164
0165 MatchingCollection match;
0166
0167 while (match.size() < partons.size()) {
0168 unsigned int partonIndex = distances[0].second / jets.size();
0169 int jetIndex = distances[0].second - jets.size() * partonIndex;
0170
0171
0172 if (useMaxDist_ && distances[0].first > maxDist_)
0173 jetIndex = -1;
0174
0175
0176 if (distances.empty())
0177 match.push_back(std::make_pair(partonIndex, -1));
0178 else
0179 match.push_back(std::make_pair(partonIndex, jetIndex));
0180
0181
0182
0183 for (unsigned int a = 0; a < distances.size(); ++a) {
0184 unsigned int pIndex = distances[a].second / jets.size();
0185 int jIndex = distances[a].second - jets.size() * pIndex;
0186 if ((pIndex == partonIndex) || (jIndex == jetIndex)) {
0187 distances.erase(distances.begin() + a, distances.begin() + a + 1);
0188 --a;
0189 }
0190 }
0191 }
0192
0193 matching.clear();
0194 matching.push_back(match);
0195 return;
0196 }
0197
0198 void JetPartonMatching::minSumDist_recursion(const unsigned int ip,
0199 std::vector<unsigned int>& jetIndices,
0200 std::vector<bool>& usedJets,
0201 std::vector<std::pair<double, MatchingCollection> >& distMatchVec) {
0202
0203 if (ip < partons.size()) {
0204 for (unsigned int ij = 0; ij < jets.size(); ++ij) {
0205 if (usedJets[ij])
0206 continue;
0207 usedJets[ij] = true;
0208 jetIndices[ip] = ij;
0209 minSumDist_recursion(ip + 1, jetIndices, usedJets, distMatchVec);
0210 usedJets[ij] = false;
0211 }
0212 return;
0213 }
0214
0215
0216 double sumDist = 0;
0217 MatchingCollection match;
0218 for (unsigned int ip = 0; ip < partons.size(); ++ip) {
0219 double dist = distance(partons[ip]->p4(), jets[jetIndices[ip]]->p4());
0220 if (useMaxDist_ && dist > maxDist_)
0221 return;
0222 sumDist += distance(partons[ip]->p4(), jets[jetIndices[ip]]->p4());
0223 match.push_back(std::make_pair(ip, jetIndices[ip]));
0224 }
0225
0226 distMatchVec.push_back(std::make_pair(sumDist, match));
0227 return;
0228 }
0229
0230 void JetPartonMatching::matchingMinSumDist() {
0231
0232
0233
0234 std::vector<std::pair<double, MatchingCollection> > distMatchVec;
0235
0236 std::vector<bool> usedJets;
0237 for (unsigned int i = 0; i < jets.size(); ++i) {
0238 usedJets.push_back(false);
0239 }
0240
0241 std::vector<unsigned int> jetIndices;
0242 jetIndices.reserve(partons.size());
0243
0244 minSumDist_recursion(0, jetIndices, usedJets, distMatchVec);
0245
0246 std::sort(distMatchVec.begin(), distMatchVec.end());
0247
0248 matching.clear();
0249
0250 if (distMatchVec.empty()) {
0251 MatchingCollection dummyMatch;
0252 for (unsigned int ip = 0; ip < partons.size(); ++ip)
0253 dummyMatch.push_back(std::make_pair(ip, -1));
0254 matching.push_back(dummyMatch);
0255 } else
0256 for (unsigned int i = 0; i < distMatchVec.size(); ++i)
0257 matching.push_back(distMatchVec[i].second);
0258
0259 return;
0260 }
0261
0262 void JetPartonMatching::matchingPtOrderedMinDist() {
0263
0264
0265
0266 std::vector<std::pair<double, unsigned int> > ptOrderedPartons;
0267
0268 for (unsigned int ip = 0; ip < partons.size(); ++ip)
0269 ptOrderedPartons.push_back(std::make_pair(partons[ip]->pt(), ip));
0270
0271 std::sort(ptOrderedPartons.begin(), ptOrderedPartons.end());
0272 std::reverse(ptOrderedPartons.begin(), ptOrderedPartons.end());
0273
0274 std::vector<unsigned int> jetIndices;
0275 for (unsigned int ij = 0; ij < jets.size(); ++ij)
0276 jetIndices.push_back(ij);
0277
0278 MatchingCollection match;
0279
0280 for (unsigned int ip = 0; ip < ptOrderedPartons.size(); ++ip) {
0281 double minDist = 999.;
0282 int ijMin = -1;
0283
0284 for (unsigned int ij = 0; ij < jetIndices.size(); ++ij) {
0285 double dist = distance(partons[ptOrderedPartons[ip].second]->p4(), jets[jetIndices[ij]]->p4());
0286 if (dist < minDist) {
0287 if (!useMaxDist_ || dist <= maxDist_) {
0288 minDist = dist;
0289 ijMin = ij;
0290 }
0291 }
0292 }
0293
0294 if (ijMin >= 0) {
0295 match.push_back(std::make_pair(ptOrderedPartons[ip].second, jetIndices[ijMin]));
0296 jetIndices.erase(jetIndices.begin() + ijMin, jetIndices.begin() + ijMin + 1);
0297 } else
0298 match.push_back(std::make_pair(ptOrderedPartons[ip].second, -1));
0299 }
0300
0301 matching.clear();
0302 matching.push_back(match);
0303 return;
0304 }
0305
0306 void JetPartonMatching::matchingUnambiguousOnly() {
0307
0308
0309 std::vector<bool> jetMatched;
0310 for (unsigned int ij = 0; ij < jets.size(); ++ij)
0311 jetMatched.push_back(false);
0312
0313 MatchingCollection match;
0314
0315 for (unsigned int ip = 0; ip < partons.size(); ++ip) {
0316 int iMatch = -1;
0317 for (unsigned int ij = 0; ij < jets.size(); ++ij) {
0318 double dist = distance(partons[ip]->p4(), jets[ij]->p4());
0319 if (dist <= maxDist_) {
0320 if (!jetMatched[ij]) {
0321 jetMatched[ij] = true;
0322 if (iMatch == -1)
0323 iMatch = ij;
0324 else
0325 iMatch = -2;
0326 } else
0327 iMatch = -2;
0328 }
0329 }
0330 match.push_back(std::make_pair(ip, iMatch));
0331 }
0332
0333 matching.clear();
0334 matching.push_back(match);
0335 return;
0336 }
0337
0338 int JetPartonMatching::getMatchForParton(const unsigned int part, const unsigned int comb) {
0339
0340
0341
0342 if (comb >= matching.size())
0343 return -9;
0344 if (part >= matching[comb].size())
0345 return -9;
0346 return (matching[comb])[part].second;
0347 }
0348
0349 std::vector<int> JetPartonMatching::getMatchesForPartons(const unsigned int comb) {
0350
0351
0352 std::vector<int> jetIndices;
0353 for (unsigned int part = 0; part < partons.size(); ++part)
0354 jetIndices.push_back(getMatchForParton(part, comb));
0355 return jetIndices;
0356 }
0357
0358 double JetPartonMatching::getDistanceForParton(const unsigned int part, const unsigned int comb) {
0359
0360 if (getMatchForParton(part, comb) < 0)
0361 return -999.;
0362 return distance(jets[getMatchForParton(part, comb)]->p4(), partons[part]->p4());
0363 }
0364
0365 double JetPartonMatching::getSumDistances(const unsigned int comb) {
0366
0367 double sumDists = 0.;
0368 for (unsigned int part = 0; part < partons.size(); ++part) {
0369 double dist = getDistanceForParton(part, comb);
0370 if (dist < 0.)
0371 return -999.;
0372 sumDists += dist;
0373 }
0374 return sumDists;
0375 }
0376
0377 void JetPartonMatching::print() {
0378
0379 edm::LogInfo log("JetPartonMatching");
0380 log << "++++++++++++++++++++++++++++++++++++++++++++++ \n";
0381 log << " algorithm : ";
0382 switch (algorithm_) {
0383 case totalMinDist:
0384 log << "totalMinDist ";
0385 break;
0386 case minSumDist:
0387 log << "minSumDist ";
0388 break;
0389 case ptOrderedMinDist:
0390 log << "ptOrderedMinDist";
0391 break;
0392 case unambiguousOnly:
0393 log << "unambiguousOnly ";
0394 break;
0395 default:
0396 log << "UNKNOWN ";
0397 }
0398 log << "\n";
0399 log << " useDeltaR : " << std::boolalpha << useDeltaR_;
0400 log << "\n";
0401 log << " useMaxDist: " << std::boolalpha << useMaxDist_;
0402 log << " maxDist: " << maxDist_ << "\n";
0403 log << " number of partons / jets: " << partons.size() << " / " << jets.size() << "\n";
0404 log << " number of available combinations: " << getNumberOfAvailableCombinations() << "\n";
0405 for (unsigned int comb = 0; comb < matching.size(); ++comb) {
0406 log << " -------------------------------------------- \n";
0407 log << " ind. of matched jets:";
0408 for (unsigned int part = 0; part < partons.size(); ++part)
0409 log << std::setw(4) << getMatchForParton(part, comb);
0410 log << "\n";
0411 log << " sumDeltaR : " << getSumDeltaR(comb) << "\n";
0412 log << " sumDeltaPt / sumDeltaE: " << getSumDeltaPt(comb) << " / " << getSumDeltaE(comb);
0413 log << "\n";
0414 }
0415 log << "++++++++++++++++++++++++++++++++++++++++++++++";
0416 }