TPS.h

Go to the documentation of this file.

#pragma once
 
#include "Core/Objects/interface/Jet.h"
#include "Core/Objects/interface/Di.h"
 
#include <array>
#include <algorithm>
#include <list>
#include <vector>
#include <numeric>
#include "Math/Vector2D.h"
#include "Math/Vector4D.h"
#include "Math/VectorUtil.h"
 
#include "Core/JetObservables/bin/common.h"
 
#include <darwin.h>
 
namespace DE = Darwin::Exceptions;
 
namespace DAS::TPS {
 
template<class Jet>
struct MultiJets {
 
    std::vector<Jet> jets;
    Di<const Jet, const Jet> leading_dijet;
 
    std::array<std::array<float, 6>, 6> pt_pair_asymm_vec;
    std::array<std::array<float, 6>, 6> pt_pair_asymm_abs;
    std::array<std::array<float, 6>, 6> Delta_eta_pair;
    std::array<std::array<float, 6>, 6> Delta_phi_pair;
    std::array<std::array<float, 6>, 6> m_inv;
 
    float S_pt = 0, S_pt_min = 0;
    std::array<float, 15> S_pt_perm{};
    std::array<float, 3> S_pt_4{}, S_pt_4_min{};
    std::array<std::array<float, 3>, 3> S_pt_4_perm;
 
    float S_phi = 0, S_phi_min = 0;
    std::array<float, 15> S_phi_perm{};
    std::array<float, 3> S_phi_4{}, S_phi_4_min{};
    std::array<std::array<float, 3>, 3> S_phi_4_perm;
 
    MultiJets(const std::vector<Jet>& inputJets, int minpt):
    jets(PhaseSpaceSelection(inputJets, minpt))
    {
        if (jets.size() > 1){
            leading_dijet = make_di(jets.at(0), jets.at(1));
        }
        if (jets.size() >= 6) FillFeatures();
    }
 
    // At least 2 jets
    explicit operator bool() const {
        return jets.size() > 1;
    }
 
    double jetswgt() const { 
        return accumulate(jets.begin(), jets.end(), 1.0, 
        [] (float w, const Jet& j) {
            return w * j.weights.front(); 
        }); 
    }
 
    float Ht () const
    {
        return accumulate(jets.begin(), jets.end(), 0.0,
            [](float sum, const auto& jet) {
                return sum + jet.CorrPt();
            });
    }
 
private:
    void FillFeatures() {
        LoopOverPairs();
        LoopOverTriplets();
        LoopOverQuartets();
    }
 
    void LoopOverPairs() {
        for (size_t i = 0; i < 6; ++i) {
            for (size_t j = 0; j < 6; ++j) {
                pt_pair_asymm_abs[i][j] = abs(jets[i].CorrPt() - jets[j].CorrPt()) / (jets[i].CorrPt() + jets[j].CorrPt());
                auto dijet = make_di(jets[i], jets[j]);
                pt_pair_asymm_vec[i][j] = dijet.CorrPt() / (jets[i].CorrPt() + jets[j].CorrPt());
                Delta_eta_pair[i][j] = dijet.DeltaEta();
                Delta_phi_pair[i][j] = dijet.DeltaPhi();
                m_inv[i][j] = dijet.CorrP4().M();
            }
        }
    }
 
    void LoopOverTriplets() {
        int i = 0;
        for (const auto& triple : jetCombinations) {
            for (auto& pr : triple) {
                auto dijet = make_di(jets[pr.first], jets[pr.second]);
                S_phi_perm[i] += pow(dijet.DeltaPhi(), 2);
                S_pt_perm[i]  += pow(dijet.CorrPt() / (jets[pr.first].CorrPt() + jets[pr.second].CorrPt()), 2);
            }
            S_phi_perm[i] = sqrt(S_phi_perm[i] / 3.0);
            S_pt_perm[i]  = sqrt(S_pt_perm[i]  / 3.0);
            i++;
        }
        S_phi = std::accumulate(S_phi_perm.begin(), S_phi_perm.end(), 0.0) / static_cast<float>(S_phi_perm.size());
        S_pt  = std::accumulate(S_pt_perm.begin(),  S_pt_perm.end(), 0.0)  / static_cast<float>(S_pt_perm.size());
        S_phi_min = *std::min_element(S_phi_perm.begin(), S_phi_perm.end());
        S_pt_min  = *std::min_element( S_pt_perm.begin(),  S_pt_perm.end());
    }
 
    void LoopOverQuartets() {
        std::array<std::array<int,4>, 3> group4 = {{
            {{0, 1, 2, 3}},
            {{0, 1, 4, 5}},
            {{2, 3, 4, 5}}
        }};
        for (size_t j = 0; j < group4.size(); ++j) {
            auto partitions = GetJetQuadPartitions(group4[j]);
            for (size_t i = 0; i < partitions.size(); ++i) {
                auto [p1, p2] = partitions[i];
                auto dijet1 = make_di(jets[p1.first], jets[p1.second]);
                auto dijet2 = make_di(jets[p2.first], jets[p2.second]);
                float a = dijet1.CorrPt() / (jets[p1.first].CorrPt() + jets[p1.second].CorrPt());
                float b = dijet2.CorrPt() / (jets[p2.first].CorrPt() + jets[p2.second].CorrPt());
                S_pt_4_perm[j][i]  = std::sqrt(0.5 * (pow(a, 2) + pow(b, 2)));
                S_phi_4_perm[j][i] = std::sqrt(0.5 * (pow(dijet1.DeltaPhi(), 2) + pow(dijet2.DeltaPhi(), 2)));
            }
            S_pt_4[j] = std::accumulate(S_pt_4_perm[j].begin(), S_pt_4_perm[j].end(), 0.0) / static_cast<float>(S_pt_4_perm[j].size());
            S_phi_4[j] = std::accumulate(S_phi_4_perm[j].begin(), S_phi_4_perm[j].end(), 0.0) / static_cast<float>(S_phi_4_perm[j].size());
            S_pt_4_min[j] = *std::min_element(S_pt_4_perm[j].begin(), S_pt_4_perm[j].end());
            S_phi_4_min[j] = *std::min_element(S_phi_4_perm[j].begin(), S_phi_4_perm[j].end());
        }
        
    }
 
    // Generates all possible partitions of 4 jets into 2 pairs
    std::array<std::pair<std::pair<int, int>, std::pair<int, int>>, 3>
    GetJetQuadPartitions(const std::array<int, 4>& index) {
        return {{
            {{index[0], index[1]}, {index[2], index[3]}},
            {{index[0], index[2]}, {index[1], index[3]}},
            {{index[0], index[3]}, {index[1], index[2]}}
        }};
    }
 
    // Generates all 15 unique pairings of 6 elements (labelled 0..5),
    // Each result is a vector of 3 std::pairs<int,int>.
    static const std::vector<std::vector<std::pair<int, int>>>& getJetPairings() {
        static const auto jetPairs = []() {
            std::vector<int> elements{0, 1, 2, 3, 4, 5};
            std::vector<std::pair<int, int>> current;
            std::vector<std::vector<std::pair<int, int>>> results;
 
            // Helper function that recursively generates all ways to partition 'elements' 
            // into pairs. Each unique partition is pushed into 'allPairs' as a list of
            // (a,b) pairs, stored in current recursion via 'current'.
            auto pairJetsIndices = [](const std::vector<int>& elements,
                                      std::vector<std::pair<int, int>>& current,
                                      std::vector<std::vector<std::pair<int, int>>>& results,
                                      auto&& self) -> void {
                // If no elements remain, we've formed all pairs for this arrangement
                if (elements.empty()) {
                    results.push_back(current);
                    return;
                }
                // Take the first element
                int first = elements.front();
                // Remaining elements
                std::vector<int> rest(elements.begin() + 1, elements.end());
 
                // Try pairing 'first' with each of the others
                for (size_t i = 0; i < rest.size(); ++i) {
                    int partner = rest[i];
                    // Add this pair to 'current'
                    current.emplace_back(first, partner);
 
                    // Build the list for the next recursion by removing 'partner'
                    std::vector<int> next;
                    next.reserve(rest.size() - 1);
                    for (size_t j = 0; j < rest.size(); ++j) {
                        if (j != i)
                            next.push_back(rest[j]);
                    }
 
                    // Recurse with the remaining elements
                    self(next, current, results, self);
                    current.pop_back();
                }
            };
            // Start the recursion
            pairJetsIndices(elements, current, results, pairJetsIndices);
            return results;
        }();
        return jetPairs;
    }
 
    // `jetCombinations` stores all possible unique ways of dividing six elements into three pairs. 
    const std::vector<std::vector<std::pair<int, int>>>& jetCombinations = getJetPairings();
 
};
 
} // end of namespace DAS::TPS