File indexing completed on 2023-03-17 11:14:55
0001 #include "MuonAnalysis/MuonAssociators/interface/MatcherByPullsAlgorithm.h"
0002
0003
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);
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 }
