JMEmatching.h

Go to the documentation of this file.

#pragma once
 
#include <iostream>
#include <queue>
#include <utility>
#include <vector>
 
#include "Core/Objects/interface/Jet.h"
 
#include "Math/VectorUtil.h"
 
namespace DAS {
 
template<typename RecContainer = std::vector<RecJet>,
         typename GenContainer = std::vector<GenJet>> class JMEmatching {
 
    static FourVector getP4 (const FourVector& p4)
    {
        return p4;
    }
 
    static FourVector getP4 (const PhysicsObject& phys_obj)
    {
        return phys_obj.p4;
    }
 
public:
 
    using RecIt = RecContainer::iterator;
    using GenIt = GenContainer::iterator;
    using Couple = std::pair<RecIt, GenIt>;
 
    static float maxDR;
 
    std::vector<RecIt> fake_its;
    std::vector<GenIt> miss_its;
    std::vector<Couple> match_its;
 
    JMEmatching (RecContainer& recJets, 
                 GenContainer& genJets) 
    {
        using namespace std;
 
        // 1) get the size of the collections
 
        size_t nRec = distance(recJets.begin(), recJets.end()),
               nGen = distance(genJets.begin(), genJets.end());
 
        // 2) declare a few utilities
 
        using Indices = pair<size_t,size_t>;
 
        auto DeltaR2 = [&recJets,&genJets](const Indices& indices) {
            const FourVector& recJet = getP4(recJets[indices.first ]),
                            & genJet = getP4(genJets[indices.second]);
            return ROOT::Math::VectorUtil::DeltaR2(recJet, genJet);
        };
 
        auto ascending_DR2 = [&DeltaR2](const Indices& l, const Indices& r) {
            return DeltaR2(l) > DeltaR2(r);
        };
 
        // 3) make pair candidates
 
        // https://en.cppreference.com/w/cpp/container/priority_queue
        priority_queue<      Indices,
                      vector<Indices>, decltype(ascending_DR2)>
                                     candidates(ascending_DR2);
        for (size_t iRec = 0; iRec < nRec; ++iRec)
        for (size_t iGen = 0; iGen < nGen; ++iGen) {
            Indices indices {iRec, iGen};
            auto DR2 = DeltaR2(indices);
            const auto maxDR2 = pow(maxDR,2);
            if (DR2 >= maxDR2) continue;
            candidates.emplace( std::move(indices) );
        }
 
        // 4) keep only the closest pairs (keep track of matched jets to avoid double matching)
 
        vector<bool> matchedRec(nRec, false),
                     matchedGen(nGen, false);
        match_its.reserve(min(nRec,nGen));
        while (!candidates.empty()) {
            auto [iRec, iGen] = candidates.top();
            candidates.pop();
 
            // skip already matched objects
            if (matchedRec.at(iRec)) continue;
            if (matchedGen.at(iGen)) continue;
            matchedRec.at(iRec) = true;
            matchedGen.at(iGen) = true;
 
            RecIt rec_it = recJets.begin();
            GenIt gen_it = genJets.begin();
            advance(rec_it, iRec);
            advance(gen_it, iGen);
 
            match_its.push_back( make_pair(rec_it, gen_it) );
        }
 
        // 5) fill the unmatched entries
 
        fake_its.reserve(nRec);
        for (int iRec = 0; iRec < nRec; ++iRec) {
            if (matchedRec.at(iRec)) continue;
            RecIt fake_it = recJets.begin();
            advance(fake_it, iRec);
            fake_its.push_back(fake_it);
        }
 
        miss_its.reserve(nGen);
        for (int iGen = 0; iGen < nGen; ++iGen) {
            if (matchedGen.at(iGen)) continue;
            GenIt miss_it = genJets.begin();
            advance(miss_it, iGen);
            miss_its.push_back(miss_it);
        }
    }
};
template<typename RecContainer, typename GenContainer>
    float JMEmatching<RecContainer, GenContainer>::maxDR = 0.2;
 
} // end of DAS namespace