MatcherByPullsAlgorithm.cc

Go to the documentation of this file.

#include "MuonAnalysis/MuonAssociators/interface/MatcherByPullsAlgorithm.h"

// user include files
#include "DataFormats/Math/interface/deltaR.h"

/*     ____                _                   _             
 *    / ___|___  _ __  ___| |_ _ __ _   _  ___| |_ ___  _ __ 
 *   | |   / _ \| '_ \/ __| __| '__| | | |/ __| __/ _ \| '__|
 *   | |__| (_) | | | \__ \ |_| |  | |_| | (__| || (_) | |   
 *    \____\___/|_| |_|___/\__|_|   \__,_|\___|\__\___/|_|   
 *                                                           
 */  
MatcherByPullsAlgorithm::MatcherByPullsAlgorithm(const edm::ParameterSet &iConfig) :
    dr2_(std::pow(iConfig.getParameter<double>("maxDeltaR"),2)),
    cut_(iConfig.getParameter<double>("maxPull")),
    diagOnly_(iConfig.getParameter<bool>("diagonalElementsOnly")),
    useVertex_(iConfig.getParameter<bool>("useVertexVariables"))
{
    std::string track = iConfig.getParameter<std::string>("track");
    if      (track == "standAloneMuon") track_ = StaTrack;
    else if (track == "combinedMuon")   track_ = GlbTrack;
    else if (track == "track")          track_ = TrkTrack;
    else throw cms::Exception("Configuration") << "MatcherByPullsAlgorithm: track '"<<track<<"' is not known\n" << 
               "Allowed values are: 'track', 'combinedMuon', 'standAloneMuon' (as per reco::RecoCandidate object)\n";
}

MatcherByPullsAlgorithm::~MatcherByPullsAlgorithm()
{
}

/*    __  __       _       _                   
 *   |  \/  | __ _| |_ ___| |__   ___ _ __ ___ 
 *   | |\/| |/ _` | __/ __| '_ \ / _ \ '__/ __|
 *   | |  | | (_| | || (__| | | |  __/ |  \__ \
 *   |_|  |_|\__,_|\__\___|_| |_|\___|_|  |___/
 *                                             
 */  
std::pair<bool,float> 
MatcherByPullsAlgorithm::match(const reco::Track &tk, const reco::Candidate &c, const AlgebraicSymMatrix55 &invCov) const {
    if (::deltaR2(tk,c) <= dr2_) {
        AlgebraicVector5 diff(tk.qoverp() - c.charge()/c.p(), 
                              tk.theta() - c.theta(), 
                              ::deltaPhi(tk.phi(), c.phi()),
                              tk.dxy(c.vertex()),
                              tk.dsz(c.vertex())); 
        double pull = ROOT::Math::Similarity(diff, invCov);
#if 0
        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;
        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;
        std::cout << "Delta: " << diff << std::endl;
        std::cout << "Sigmas: ";
        for (size_t i = 0; i < 5; ++i) {
            if (invCov(i,i) == 0) std::cout << "---\t";
            else std::cout << std::sqrt(1.0/invCov(i,i)) << "\t";
        }
        std::cout << std::endl;
        std::cout << "Items: ";
        for (size_t i = 0; i < 5; ++i) {
            if (invCov(i,i) == 0) std::cout << "---\t";
            else std::cout << diff(i)*std::sqrt(invCov(i,i)) << "\t";
        }
        std::cout << std::endl;
        std::cout << "Pull: "  << pull << std::endl;
#endif
        return std::pair<bool,float>(pull < cut_, pull);
    }
    return std::pair<bool,float>(false,9e9);
}



std::pair<int,float>
MatcherByPullsAlgorithm::match(const reco::RecoCandidate &src, 
                      const std::vector<reco::GenParticle> &cands,
                      const std::vector<uint8_t>       &good) const 
{
    const reco::Track * tk = track(src);
    return (tk == nullptr ? 
            std::pair<int,float>(-1,9e9) : 
            match(*tk, cands, good));
}

std::pair<int,float>
MatcherByPullsAlgorithm::match(const reco::Track &tk, 
                      const std::vector<reco::GenParticle> &cands,
                      const std::vector<uint8_t>       &good) const 
{
    std::pair<int,float> best(-1,9e9);

    AlgebraicSymMatrix55 invCov; fillInvCov(tk, invCov);
    for (int i = 0, n = cands.size(); i < n; ++i) {
        if (!good[i]) continue;
        std::pair<bool,float> m = match(tk, cands[i], invCov);
        if (m.first && (m.second < best.second)) {
            best.first  = i;
            best.second = m.second;
        } 
    }
    return best;
}

void
MatcherByPullsAlgorithm::matchMany(const reco::RecoCandidate &src,
                const std::vector<reco::GenParticle> &cands,
                const std::vector<uint8_t>           &good,
                std::vector<std::pair<double,int> >  &matchesToFill) const
{
    const reco::Track * tk = track(src);
    if (tk != nullptr) matchMany(*tk, cands, good, matchesToFill);
}

void
MatcherByPullsAlgorithm::matchMany(const reco::Track &tk,
                const std::vector<reco::GenParticle> &cands,
                const std::vector<uint8_t>           &good,
                std::vector<std::pair<double,int> >  &matchesToFill) const
{

    AlgebraicSymMatrix55 invCov; fillInvCov(tk, invCov);
    for (int i = 0, n = cands.size(); i < n; ++i) {
        if (!good[i]) continue;
        std::pair<bool,double> m = match(tk, cands[i], invCov);
        if (m.first) matchesToFill.push_back(std::make_pair(m.second,i));
    }
    std::sort(matchesToFill.begin(),matchesToFill.end());
}

/*    _   _ _   _ _ _ _   _           
 *   | | | | |_(_) (_) |_(_) ___  ___ 
 *   | | | | __| | | | __| |/ _ \/ __|
 *   | |_| | |_| | | | |_| |  __/\__ \
 *    \___/ \__|_|_|_|\__|_|\___||___/
 *                                    
 */  
const reco::Track * 
MatcherByPullsAlgorithm::track(const reco::RecoCandidate &muon) const {
    switch (track_) {
        case StaTrack : return muon.standAloneMuon().isNonnull() ? muon.standAloneMuon().get() : nullptr;
        case GlbTrack : return muon.combinedMuon().isNonnull()   ? muon.combinedMuon().get()   : nullptr;
        case TrkTrack : return muon.track().isNonnull()          ? muon.track().get()          : nullptr;
    }
    assert(false);
}


void
MatcherByPullsAlgorithm::fillInvCov(const reco::Track &tk, AlgebraicSymMatrix55 &invCov) const 
{
    if (useVertex_) {
        invCov = tk.covariance();
        if (diagOnly_) {
            for (size_t i = 0; i < 5; ++i) { for (size_t j = i+1; j < 5; ++j) { 
                invCov(i,j) = 0;
            } }
        }
        invCov.Invert();
    } else {
        AlgebraicSymMatrix33 momCov = tk.covariance().Sub<AlgebraicSymMatrix33>(0,0); // get 3x3 matrix
        if (diagOnly_) { momCov(0,1) = 0; momCov(0,2) = 0; momCov(1,2) = 0; }
        momCov.Invert();
        invCov.Place_at(momCov,0,0);
    }
}