Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 13:31:43

0001 #include "MuonAnalysis/MuonAssociators/interface/MatcherByPullsAlgorithm.h"
0002 
0003 // user include files
0004 #include "DataFormats/Math/interface/deltaR.h"
0005 
0006 /*     ____                _                   _             
0007  *    / ___|___  _ __  ___| |_ _ __ _   _  ___| |_ ___  _ __ 
0008  *   | |   / _ \| '_ \/ __| __| '__| | | |/ __| __/ _ \| '__|
0009  *   | |__| (_) | | | \__ \ |_| |  | |_| | (__| || (_) | |   
0010  *    \____\___/|_| |_|___/\__|_|   \__,_|\___|\__\___/|_|   
0011  *                                                           
0012  */
0013 MatcherByPullsAlgorithm::MatcherByPullsAlgorithm(const edm::ParameterSet &iConfig)
0014     : dr2_(std::pow(iConfig.getParameter<double>("maxDeltaR"), 2)),
0015       cut_(iConfig.getParameter<double>("maxPull")),
0016       diagOnly_(iConfig.getParameter<bool>("diagonalElementsOnly")),
0017       useVertex_(iConfig.getParameter<bool>("useVertexVariables")) {
0018   std::string track = iConfig.getParameter<std::string>("track");
0019   if (track == "standAloneMuon")
0020     track_ = StaTrack;
0021   else if (track == "combinedMuon")
0022     track_ = GlbTrack;
0023   else if (track == "track")
0024     track_ = TrkTrack;
0025   else
0026     throw cms::Exception("Configuration")
0027         << "MatcherByPullsAlgorithm: track '" << track << "' is not known\n"
0028         << "Allowed values are: 'track', 'combinedMuon', 'standAloneMuon' (as per reco::RecoCandidate object)\n";
0029 }
0030 
0031 MatcherByPullsAlgorithm::~MatcherByPullsAlgorithm() {}
0032 
0033 /*    __  __       _       _                   
0034  *   |  \/  | __ _| |_ ___| |__   ___ _ __ ___ 
0035  *   | |\/| |/ _` | __/ __| '_ \ / _ \ '__/ __|
0036  *   | |  | | (_| | || (__| | | |  __/ |  \__ \
0037  *   |_|  |_|\__,_|\__\___|_| |_|\___|_|  |___/
0038  *                                             
0039  */
0040 std::pair<bool, float> MatcherByPullsAlgorithm::match(const reco::Track &tk,
0041                                                       const reco::Candidate &c,
0042                                                       const AlgebraicSymMatrix55 &invCov) const {
0043   if (::deltaR2(tk, c) <= dr2_) {
0044     AlgebraicVector5 diff(tk.qoverp() - c.charge() / c.p(),
0045                           tk.theta() - c.theta(),
0046                           ::deltaPhi(tk.phi(), c.phi()),
0047                           tk.dxy(c.vertex()),
0048                           tk.dsz(c.vertex()));
0049     double pull = ROOT::Math::Similarity(diff, invCov);
0050 #if 0
0051         std::cout << "Tk charge/pt/eta/phi/vx/vy/vz " << tk.charge() << "\t" << tk.pt() << "\t" << tk.eta() << "\t" << tk.phi() << "\t" << tk.vx() << "\t" << tk.vy() << "\t" << tk.vz() << std::endl;
0052         std::cout << "MC charge/pt/eta/phi/vx/vy/vz " << c.charge() << "\t" << c.pt() << "\t" << c.eta() << "\t" << c.phi() << "\t" << c.vx() << "\t" << c.vy() << "\t" << c.vz() << std::endl;
0053         std::cout << "Delta: " << diff << std::endl;
0054         std::cout << "Sigmas: ";
0055         for (size_t i = 0; i < 5; ++i) {
0056             if (invCov(i,i) == 0) std::cout << "---\t";
0057             else std::cout << std::sqrt(1.0/invCov(i,i)) << "\t";
0058         }
0059         std::cout << std::endl;
0060         std::cout << "Items: ";
0061         for (size_t i = 0; i < 5; ++i) {
0062             if (invCov(i,i) == 0) std::cout << "---\t";
0063             else std::cout << diff(i)*std::sqrt(invCov(i,i)) << "\t";
0064         }
0065         std::cout << std::endl;
0066         std::cout << "Pull: "  << pull << std::endl;
0067 #endif
0068     return std::pair<bool, float>(pull < cut_, pull);
0069   }
0070   return std::pair<bool, float>(false, 9e9);
0071 }
0072 
0073 std::pair<int, float> MatcherByPullsAlgorithm::match(const reco::RecoCandidate &src,
0074                                                      const std::vector<reco::GenParticle> &cands,
0075                                                      const std::vector<uint8_t> &good) const {
0076   const reco::Track *tk = track(src);
0077   return (tk == nullptr ? std::pair<int, float>(-1, 9e9) : match(*tk, cands, good));
0078 }
0079 
0080 std::pair<int, float> MatcherByPullsAlgorithm::match(const reco::Track &tk,
0081                                                      const std::vector<reco::GenParticle> &cands,
0082                                                      const std::vector<uint8_t> &good) const {
0083   std::pair<int, float> best(-1, 9e9);
0084 
0085   AlgebraicSymMatrix55 invCov;
0086   fillInvCov(tk, invCov);
0087   for (int i = 0, n = cands.size(); i < n; ++i) {
0088     if (!good[i])
0089       continue;
0090     std::pair<bool, float> m = match(tk, cands[i], invCov);
0091     if (m.first && (m.second < best.second)) {
0092       best.first = i;
0093       best.second = m.second;
0094     }
0095   }
0096   return best;
0097 }
0098 
0099 void MatcherByPullsAlgorithm::matchMany(const reco::RecoCandidate &src,
0100                                         const std::vector<reco::GenParticle> &cands,
0101                                         const std::vector<uint8_t> &good,
0102                                         std::vector<std::pair<double, int> > &matchesToFill) const {
0103   const reco::Track *tk = track(src);
0104   if (tk != nullptr)
0105     matchMany(*tk, cands, good, matchesToFill);
0106 }
0107 
0108 void MatcherByPullsAlgorithm::matchMany(const reco::Track &tk,
0109                                         const std::vector<reco::GenParticle> &cands,
0110                                         const std::vector<uint8_t> &good,
0111                                         std::vector<std::pair<double, int> > &matchesToFill) const {
0112   AlgebraicSymMatrix55 invCov;
0113   fillInvCov(tk, invCov);
0114   for (int i = 0, n = cands.size(); i < n; ++i) {
0115     if (!good[i])
0116       continue;
0117     std::pair<bool, double> m = match(tk, cands[i], invCov);
0118     if (m.first)
0119       matchesToFill.push_back(std::make_pair(m.second, i));
0120   }
0121   std::sort(matchesToFill.begin(), matchesToFill.end());
0122 }
0123 
0124 /*    _   _ _   _ _ _ _   _           
0125  *   | | | | |_(_) (_) |_(_) ___  ___ 
0126  *   | | | | __| | | | __| |/ _ \/ __|
0127  *   | |_| | |_| | | | |_| |  __/\__ \
0128  *    \___/ \__|_|_|_|\__|_|\___||___/
0129  *                                    
0130  */
0131 const reco::Track *MatcherByPullsAlgorithm::track(const reco::RecoCandidate &muon) const {
0132   switch (track_) {
0133     case StaTrack:
0134       return muon.standAloneMuon().isNonnull() ? muon.standAloneMuon().get() : nullptr;
0135     case GlbTrack:
0136       return muon.combinedMuon().isNonnull() ? muon.combinedMuon().get() : nullptr;
0137     case TrkTrack:
0138       return muon.track().isNonnull() ? muon.track().get() : nullptr;
0139   }
0140   assert(false);
0141 }
0142 
0143 void MatcherByPullsAlgorithm::fillInvCov(const reco::Track &tk, AlgebraicSymMatrix55 &invCov) const {
0144   if (useVertex_) {
0145     invCov = tk.covariance();
0146     if (diagOnly_) {
0147       for (size_t i = 0; i < 5; ++i) {
0148         for (size_t j = i + 1; j < 5; ++j) {
0149           invCov(i, j) = 0;
0150         }
0151       }
0152     }
0153     invCov.Invert();
0154   } else {
0155     AlgebraicSymMatrix33 momCov = tk.covariance().Sub<AlgebraicSymMatrix33>(0, 0);  // get 3x3 matrix
0156     if (diagOnly_) {
0157       momCov(0, 1) = 0;
0158       momCov(0, 2) = 0;
0159       momCov(1, 2) = 0;
0160     }
0161     momCov.Invert();
0162     invCov.Place_at(momCov, 0, 0);
0163   }
0164 }