TwoBody.hh

Go to the documentation of this file.

#pragma once

#include <boost/range/algorithm/unique.hpp>
#include <boost/range/algorithm_ext/erase.hpp>

#include "boca/fastjet/ClusterSequence.hh"
#include "boca/generic/Vector.hh"
#include "boca/multiplets/Multiplet.hh"
#include "boca/fastjet/Sort.hh"
#include "boca/Settings.hh"

namespace boca
{

template <typename Multiplet_1_, typename Multiplet_2_>
class TwoBody : public Multiplet
{

public:

    TwoBody() {};

    TwoBody(Multiplet_1_ const &multiplet_1, Multiplet_2_ const &multiplet_2)
    {
        multiplet_1_ = multiplet_1;
        multiplet_2_ = multiplet_2;
        if (multiplet_1.Bdt() != InitialValue() && multiplet_2.Bdt() != InitialValue()) SetBdt(multiplet_1.Bdt(), multiplet_2.Bdt());
        else if (multiplet_1.Bdt() != InitialValue()) SetBdt(multiplet_1.Bdt());
        else if (multiplet_2.Bdt() != InitialValue()) SetBdt(multiplet_2.Bdt());
    }

    template<typename Multiplet_3_, typename Multiplet_4_>
    void Enforce(TwoBody<Multiplet_3_, Multiplet_4_> const &multiplet)
    {
        multiplet_1_.Enforce(multiplet.Multiplet1());
        multiplet_2_.Enforce(multiplet.Multiplet2());
        SetBdt(multiplet.Bdt());
    }

    void Enforce(boca::Jet const &jet)
    {
        Enforce(jet, jet.Bdt());
    }

    void Enforce(std::vector<boca::Jet> const &jets)
    {
        if (jets.empty()) return;
        if (jets.size() == 1) return EnforceJet(jets.front(), jets.front().Bdt());
        multiplet_1_.Enforce(jets.at(0));
        multiplet_2_.Enforce(jets.at(1));
        SetBdt(jets.at(0).Bdt(), jets.at(1).Bdt());
        if (jets.size() > 2) std::cout << "to many jets to enforce a multiplet" << '\n';
    }

    void Enforce(boca::Jet const &jet, double bdt)
    {
        if (jet.Constituents().size() < 2) return EnforceJet(jet, bdt);
        ClusterSequence cluster_sequence(jet.Constituents(), Settings::SubJetDefinition());
        auto jets = cluster_sequence.ExclusiveJetsUpTo(2);
        if (jets.size() != 2) std::cout << "not the sub-jet number we expected: " << jets.size();
        multiplet_1_.Enforce(jets.at(0), bdt);
        multiplet_2_.Enforce(jets.at(1), bdt);
        SetBdt(bdt);
    }

    void EnforceJet(boca::Jet jet, double bdt)
    {
        jet.Info().SetSubStructure(false);
        multiplet_1_.Enforce(jet / 2, bdt);
        multiplet_2_.Enforce(jet / 2, bdt);
        SetBdt(bdt);
    }

    void Set(Multiplet_1_ const &multiplet_1, Multiplet_2_ const &multiplet_2)
    {
        multiplet_1_ = multiplet_1;
        multiplet_2_ = multiplet_2;
    }

    boca::EventShapes EventShapes() const
    {
        return event_shapes_.Get([this]() {
            return boca::EventShapes(Jets());
        });
    }

    std::vector<boca::Jet> Jets() const override
    {
        return Combine(Multiplet1().Jets(), Multiplet2().Jets());
    }

    std::vector<LorentzVector<Momentum>> LorentzVectors() const override
    {
        return DeltaR() > Settings::JetConeSize() ? Combine(Multiplet1().LorentzVectors(), Multiplet2().LorentzVectors()) : std::vector<LorentzVector<Momentum>> {Jet().LorentzVector()};
    }

    Multiplet_1_ &Multiplet1()
    {
        return multiplet_1_;
    }

    Multiplet_1_ const &Multiplet1() const
    {
        return multiplet_1_;
    }

    Multiplet_2_ &Multiplet2()
    {
        return multiplet_2_;
    }

    Multiplet_2_ const &Multiplet2() const
    {
        return multiplet_2_;
    }

    template <typename Multiplet_3_>
    bool Overlap(Multiplet_3_ const &multiplet) const
    {
        return multiplet.Overlap(multiplet_1_) || multiplet.Overlap(multiplet_2_);
    }

    bool Overlap(boca::Singlet const &singlet) const
    {
        return multiplet_1_.Overlap(singlet) || multiplet_2_.Overlap(singlet);
    }

    bool Overlap(boca::Jet const &jet) const
    {
        return multiplet_1_.Overlap(jet) || multiplet_2_.Overlap(jet);
    }

    bool Overlap() const
    {
        return multiplet_1_.Overlap(multiplet_2_);
    }

    Momentum DeltaPt() const
    {
        return Multiplet2().Pt() - Multiplet1().Pt();
    }

    Momentum Ht() const
    {
        return Multiplet1().Ht() + Multiplet2().Ht();
    }

    boca::Mass DeltaM() const
    {
        return Multiplet2().Mass() - Multiplet1().Mass();
    }

    Momentum DeltaHt() const
    {
        return Multiplet2().Ht() - Multiplet1().Ht();
    }

    boca::Mass MassDifferenceTo(Id id) const
    {
        return boost::units::abs(Mass() - MassOf(id));
    }

    Angle DeltaRap() const
    {
        return Multiplet1().DeltaRapTo(Multiplet2());
    }

    Angle DeltaPhi() const
    {
        return Multiplet1().DeltaPhiTo(Multiplet2());
    }

    Angle DeltaR() const
    {
        return Multiplet1().DeltaRTo(Multiplet2());
    }

    Line2<Angle> AngleLine() const
    {
        return {Multiplet2().Angles(), Multiplet1().AnglesMinTo(Multiplet2())};
    }

    template<typename Multiplet_>
    Angle PullTo(Multiplet_ const &multiplet) const
    {
        return ConstituentJet().PullAngle(DeltaTo(multiplet));
    }

    Angle Pull12() const
    {
        return Multiplet1().PullTo(Multiplet2());
    }

    Angle Pull21() const
    {
        return Multiplet2().PullTo(Multiplet1());
    }

    double Dipolarity() const
    {
        if (Pt() <= Settings::MinCellPt() || DeltaR() <= Settings::MinCellResolution()) return 0;
        auto sum = ConstituentJet().DipolaritySum(AngleLine());
        auto dipolarity = sum / Pt() / sqr(DeltaR());
        if(dipolarity > 100) std::cerr << "dipol: " << dipolarity << " sum: " <<  sum <<  " mom: " << Pt() << " deltar: " << DeltaR() <<  '\n';
        return dipolarity < mU ?  -6 : std::log10(dipolarity);
    }

    boca::Jet HarderComponent() const {
        return Multiplet1().Pt() > Multiplet2().Pt() ? Multiplet1().Jet() : Multiplet2().Jet();
    }

    boca::Jet SofterComponent() const {
        return Multiplet1().Pt() < Multiplet2().Pt() ? Multiplet1().Jet() : Multiplet2().Jet();
    }

    boca::Jet ComponentWithBetterMass(Id id) const {
        return Multiplet1().MassDifferenceTo(id) < Multiplet2().MassDifferenceTo(id) ? Multiplet1().Jet() : Multiplet2().Jet();
    }

    boca::Jet ComponentWithWorseMass(Id id) const {
        return Multiplet1().MassDifferenceTo(id) > Multiplet2().MassDifferenceTo(id) ? Multiplet1().Jet() : Multiplet2().Jet();
    }

    double Rho() const
    {
        return Pt() > Settings::MinCellPt() && DeltaR() > Settings::MinCellResolution() ? static_cast<double>(Mass() / Pt() / DeltaR() * 2_rad) : 0;
    }

    int Charge() const
    {
        return Multiplet1().Charge() + Multiplet2().Charge();
    }

    friend std::ostream & operator<<(std::ostream & stream, TwoBody const& two_body){
        stream << two_body.Name() << '\n' << two_body.multiplet_1_ << '\n' << two_body.multiplet_2_;
        return stream;
    }

protected:

    std::string Name() const override {
        return "Two-Body";
    }

    void SetMultiplet1(Multiplet_1_ const &multiplet_1)
    {
        multiplet_1_ = multiplet_1;
    }

    void SetMultiplet2(Multiplet_2_ const &multiplet_2)
    {
        multiplet_2_ = multiplet_2;
    }

private:

    boca::Jet GetJet() const override
    {
        return Join(Multiplet1(), Multiplet2());
    }

    Singlet GetConstituentJet() const override
    {
        return JoinConstituents(Multiplet1(), Multiplet2());
    }

    Multiplet_1_ multiplet_1_;

    Multiplet_2_ multiplet_2_;

};

template<typename Multiplet_1_, typename Multiplet_2_>
Jet Join(Multiplet_1_ const &multiplet_1, Multiplet_2_ const &multiplet_2)
{
    return Join(multiplet_1.Jet(), multiplet_2.Jet());
}

template<typename Multiplet_1_, typename Multiplet_2_>
boca::Singlet JoinConstituents(Multiplet_1_ const &multiplet_1, Multiplet_2_ const &multiplet_2)
{
    auto constituents = SortedByPt(Combine(multiplet_1.Constituents(), multiplet_2.Constituents()));
    boost::erase(constituents, boost::unique<boost::return_found_end>(constituents));
    return Join(constituents);
}

template<typename Multiplet_>
Jet Join(Jet const &jet, Multiplet_ const &multiplet)
{
    return Join(jet, multiplet.Jet());
}

template<typename Multiplet_>
Jet Join(Multiplet_ const &multiplet, Jet const &jet)
{
    return Join(jet, multiplet.Jet());
}

}