Back to home page

Project CMSSW displayed by LXR
 
 

    

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   // use maximal distance between objects
0059   // in case of unambiguousOnly algorithmm
0060   if (algorithm_ == unambiguousOnly)
0061     useMaxDist_ = true;
0062 
0063   // check if there are empty partons in
0064   // the vector, which happpens if the
0065   // event is not ttbar or the decay is
0066   // not as expected
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   // switch algorithm, default is to match
0076   // on the minimal sum of the distance
0077   // (if jets or a parton is empty fill match with blanks)
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   // calculate the distance between two lorentz vectors
0143   // using DeltaR(eta, phi) or normal space angle(theta, phi)
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   // match parton to jet with shortest distance
0151   // starting with the shortest distance available
0152   // apply some outlier rejection if desired
0153 
0154   // prepare vector of pairs with distances between
0155   // all partons to all jets in the input vectors
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     // use primitive outlier rejection if desired
0172     if (useMaxDist_ && distances[0].first > maxDist_)
0173       jetIndex = -1;
0174 
0175     // prevent underflow in case of too few jets
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     // remove all values for the matched parton
0182     // and the matched jet
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   // build up jet combinations recursively
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   // calculate sumDist for each completed combination
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;  // outlier rejection
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   // match partons to jets with minimal sum of
0232   // the distances between all partons and jets
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   // match partons to jets with minimal sum of
0264   // the distances between all partons and jets
0265   // order partons in pt first
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   // match partons to jets, only accept event
0308   // if there are no ambiguities
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]) {  // new match for jet
0321           jetMatched[ij] = true;
0322           if (iMatch == -1)  // new match for parton and jet
0323             iMatch = ij;
0324           else  // parton already matched: ambiguity!
0325             iMatch = -2;
0326         } else  // jet already matched: ambiguity!
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   // return index of the matched jet for a given parton
0340   // (if arguments for parton index and combinatoric index
0341   // are in the valid range)
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   // return a vector with the indices of the matched jets
0351   // (ordered according to the vector of partons)
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   // get the distance between parton and its best matched jet
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   // get sum of distances between partons and matched jets
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   //report using MessageLogger
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 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!