00001 #include "MuonAnalysis/MuonAssociators/interface/MatcherByPullsAlgorithm.h"
00002
00003
00004 #include "DataFormats/Math/interface/deltaR.h"
00005
00006
00007
00008
00009
00010
00011
00012
00013 MatcherByPullsAlgorithm::MatcherByPullsAlgorithm(const edm::ParameterSet &iConfig) :
00014 dr2_(std::pow(iConfig.getParameter<double>("maxDeltaR"),2)),
00015 cut_(iConfig.getParameter<double>("maxPull")),
00016 diagOnly_(iConfig.getParameter<bool>("diagonalElementsOnly")),
00017 useVertex_(iConfig.getParameter<bool>("useVertexVariables"))
00018 {
00019 std::string track = iConfig.getParameter<std::string>("track");
00020 if (track == "standAloneMuon") track_ = StaTrack;
00021 else if (track == "combinedMuon") track_ = GlbTrack;
00022 else if (track == "track") track_ = TrkTrack;
00023 else throw cms::Exception("Configuration") << "MatcherByPullsAlgorithm: track '"<<track<<"' is not known\n" <<
00024 "Allowed values are: 'track', 'combinedMuon', 'standAloneMuon' (as per reco::RecoCandidate object)\n";
00025 }
00026
00027 MatcherByPullsAlgorithm::~MatcherByPullsAlgorithm()
00028 {
00029 }
00030
00031
00032
00033
00034
00035
00036
00037
00038 std::pair<bool,float>
00039 MatcherByPullsAlgorithm::match(const reco::Track &tk, const reco::Candidate &c, const AlgebraicSymMatrix55 &invCov) const {
00040 if (::deltaR2(tk,c) <= dr2_) {
00041 AlgebraicVector5 diff(tk.qoverp() - c.charge()/c.p(),
00042 tk.theta() - c.theta(),
00043 ::deltaPhi(tk.phi(), c.phi()),
00044 tk.dxy(c.vertex()),
00045 tk.dsz(c.vertex()));
00046 double pull = ROOT::Math::Similarity(diff, invCov);
00047 #if 0
00048 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;
00049 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;
00050 std::cout << "Delta: " << diff << std::endl;
00051 std::cout << "Sigmas: ";
00052 for (size_t i = 0; i < 5; ++i) {
00053 if (invCov(i,i) == 0) std::cout << "---\t";
00054 else std::cout << std::sqrt(1.0/invCov(i,i)) << "\t";
00055 }
00056 std::cout << std::endl;
00057 std::cout << "Items: ";
00058 for (size_t i = 0; i < 5; ++i) {
00059 if (invCov(i,i) == 0) std::cout << "---\t";
00060 else std::cout << diff(i)*std::sqrt(invCov(i,i)) << "\t";
00061 }
00062 std::cout << std::endl;
00063 std::cout << "Pull: " << pull << std::endl;
00064 #endif
00065 return std::pair<bool,float>(pull < cut_, pull);
00066 }
00067 return std::pair<bool,float>(false,9e9);
00068 }
00069
00070
00071
00072 std::pair<int,float>
00073 MatcherByPullsAlgorithm::match(const reco::RecoCandidate &src,
00074 const std::vector<reco::GenParticle> &cands,
00075 const std::vector<uint8_t> &good) const
00076 {
00077 const reco::Track * tk = track(src);
00078 return (tk == 0 ?
00079 std::pair<int,float>(-1,9e9) :
00080 match(*tk, cands, good));
00081 }
00082
00083 std::pair<int,float>
00084 MatcherByPullsAlgorithm::match(const reco::Track &tk,
00085 const std::vector<reco::GenParticle> &cands,
00086 const std::vector<uint8_t> &good) const
00087 {
00088 std::pair<int,float> best(-1,9e9);
00089
00090 AlgebraicSymMatrix55 invCov; fillInvCov(tk, invCov);
00091 for (int i = 0, n = cands.size(); i < n; ++i) {
00092 if (!good[i]) continue;
00093 std::pair<bool,float> m = match(tk, cands[i], invCov);
00094 if (m.first && (m.second < best.second)) {
00095 best.first = i;
00096 best.second = m.second;
00097 }
00098 }
00099 return best;
00100 }
00101
00102 void
00103 MatcherByPullsAlgorithm::matchMany(const reco::RecoCandidate &src,
00104 const std::vector<reco::GenParticle> &cands,
00105 const std::vector<uint8_t> &good,
00106 std::vector<std::pair<double,int> > &matchesToFill) const
00107 {
00108 const reco::Track * tk = track(src);
00109 if (tk != 0) matchMany(*tk, cands, good, matchesToFill);
00110 }
00111
00112 void
00113 MatcherByPullsAlgorithm::matchMany(const reco::Track &tk,
00114 const std::vector<reco::GenParticle> &cands,
00115 const std::vector<uint8_t> &good,
00116 std::vector<std::pair<double,int> > &matchesToFill) const
00117 {
00118
00119 AlgebraicSymMatrix55 invCov; fillInvCov(tk, invCov);
00120 for (int i = 0, n = cands.size(); i < n; ++i) {
00121 if (!good[i]) continue;
00122 std::pair<bool,double> m = match(tk, cands[i], invCov);
00123 if (m.first) matchesToFill.push_back(std::make_pair(m.second,i));
00124 }
00125 std::sort(matchesToFill.begin(),matchesToFill.end());
00126 }
00127
00128
00129
00130
00131
00132
00133
00134
00135 const reco::Track *
00136 MatcherByPullsAlgorithm::track(const reco::RecoCandidate &muon) const {
00137 switch (track_) {
00138 case StaTrack : return muon.standAloneMuon().isNonnull() ? muon.standAloneMuon().get() : 0;
00139 case GlbTrack : return muon.combinedMuon().isNonnull() ? muon.combinedMuon().get() : 0;
00140 case TrkTrack : return muon.track().isNonnull() ? muon.track().get() : 0;
00141 }
00142 assert(false);
00143 }
00144
00145
00146 void
00147 MatcherByPullsAlgorithm::fillInvCov(const reco::Track &tk, AlgebraicSymMatrix55 &invCov) const
00148 {
00149 if (useVertex_) {
00150 invCov = tk.covariance();
00151 if (diagOnly_) {
00152 for (size_t i = 0; i < 5; ++i) { for (size_t j = i+1; j < 5; ++j) {
00153 invCov(i,j) = 0;
00154 } }
00155 }
00156 invCov.Invert();
00157 } else {
00158 AlgebraicSymMatrix33 momCov = tk.covariance().Sub<AlgebraicSymMatrix33>(0,0);
00159 if (diagOnly_) { momCov(0,1) = 0; momCov(0,2) = 0; momCov(1,2) = 0; }
00160 momCov.Invert();
00161 invCov.Place_at(momCov,0,0);
00162 }
00163 }
1.7.3