TPS.h

Go to the documentation of this file.

#pragma once
 
#include "Core/Objects/interface/Jet.h"
#include "Core/Objects/interface/Di.h"
 
#include <algorithm>
#include <list>
#include <vector>
#include <numeric>
#include "Math/Vector2D.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> dijet;
 
    MultiJets(const std::vector<Jet>& inputJets, int minpt):
    jets(PhaseSpaceSelection(inputJets, minpt))
    {
        if (jets.size() > 1){
            dijet = jets.at(0) + jets.at(1);
        }
    }
 
    // At least 2 jets
    explicit operator bool() const {
        return jets.size() > 1;
    }
 
    double jetswgt() const { 
        return accumulate(jets.begin(), jets.end(), 1.0, 
        [] (double w, const Jet& j) {
            return w * j.weights.front(); 
        }); 
    }
 
    double Ht () const
    {
        return accumulate(jets.begin(), jets.end(), 0.0,
            [](double sum, const auto& jet) {
                return sum + jet.CorrPt();
            });
    }
 
    float SpT() const
    {
        using namespace std;
        if (jets.size() < 6)
            BOOST_THROW_EXCEPTION( invalid_argument("insufficient jetMultiplicity for SpT calculation") );
        
        float sumS = 0.0f;
        for (const auto& triple : jetCombinations) {
            float S2 = 0.0f;
            // triple.size() == 3, each triple[i] is a pair<int,int>
            for (auto& pr : triple) {
                FourVector jet1 = jets[pr.first ].CorrP4();
                FourVector jet2 = jets[pr.second].CorrP4();
                ROOT::Math::XYVector pTsum(jet1.Px() + jet2.Px(), jet1.Py() + jet2.Py());
                S2 += pow( pTsum.R() / (jets[pr.first].CorrPt() + jets[pr.second].CorrPt()), 2 );
            }
            sumS += sqrt(S2 / 3.0f);
        }
        return sumS / static_cast<float>(jetCombinations.size());
    }
    
    float Sphi() const
    {
        using namespace std;
        if (jets.size() < 6)
            BOOST_THROW_EXCEPTION( invalid_argument("insufficient jetMultiplicity for Sphi calculation") );
 
        float sumS = 0.0f;
        for (const auto& triple : jetCombinations) {
            float S2 = 0.0f;
            for (auto& pr : triple) {
                FourVector jet1 = jets[pr.first ].CorrP4();
                FourVector jet2 = jets[pr.second].CorrP4();
                S2 += pow(DeltaPhi(jet1, jet2), 2);
            }
            sumS += sqrt(S2 / 3.0f);
        }
        return sumS / static_cast<float>(jetCombinations.size());
    }
private:
    // 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