kinematics.h

Go to the documentation of this file.

// This file is part of MARTY.
// 
// MARTY is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
// 
// MARTY is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU General Public License for more details.
// 
// You should have received a copy of the GNU General Public License
// along with MARTY. If not, see <https://www.gnu.org/licenses/>.

#pragma once

#include <csl.h>
#include <numeric>
#include "insertion.h"

namespace mty::wick {
    class Graph;
}

namespace mty {

    class Kinematics {

    public:


        Kinematics() = default;

        explicit
        Kinematics(std::vector<mty::Insertion> const &t_insertions);

        Kinematics(
                std::vector<mty::Insertion> const &t_insertions,
                std::vector<csl::Tensor>    const &t_momenta
                );

        Kinematics(
                std::vector<mty::Insertion> const &t_insertions,
                std::vector<size_t>         const &indices
                );

        size_t size() const {
            return insertions.size();
        }

        mty::Insertion const &insertion(size_t i) const {
            return insertions[i];
        }

        csl::Tensor const &momentum(size_t i) const {
            return momenta[i];
        }

        csl::Expr squaredMomentum(size_t i, size_t j) const {
            return squaredMomenta[squaredMomentumIndex(i, j)];
        }

        std::vector<mty::Insertion> const &getInsertions() const {
            return insertions;
        }

        std::vector<csl::Tensor> const &getMomenta() const {
            return momenta;
        }

        std::vector<csl::Tensor> getOrderedMomenta() const;

        std::vector<csl::Expr> const &getSquaredMomenta() const {
            return squaredMomenta;
        }

        std::vector<size_t> const &getIndices() const {
            return indices;
        }

        csl::Expr getDegreesOfFreedomFactor() const;

        Kinematics subset(std::vector<size_t> pos) const;

        void applyPermutation(std::vector<size_t> const &pos);

        Kinematics applyIndices(
                std::vector<size_t> const &indices
                ) const;

        void sortFromIndices();

        Kinematics alignedWith(Kinematics const &other) const;

        static void replace(
                csl::Expr        &expr,
                Kinematics const &k1,
                Kinematics const &k2
                );

        static Kinematics merge(
                Kinematics const &k1,
                Kinematics const &k2
                );

        static std::vector<size_t> defaultIndices(size_t N) {
            std::vector<size_t> indices(N);
            std::iota(begin(indices), end(indices), 1);
            return indices;
        }

        friend std::ostream &operator<<(
                std::ostream     &out,
                Kinematics const &kin
                );

    private:

        void initMomentaSquared(std::vector<size_t> const &indices);

        void addContraction(
                csl::Tensor     &p1,
                csl::Tensor     &p2, 
                csl::Expr const &res
                );
        
        void setSquaredMomentum(
                size_t           i, 
                size_t           j,
                csl::Expr const &res);

        size_t squaredMomentumIndex(size_t i, size_t j) const {
            return i * momenta.size() + j;
        }

    private:

        std::vector<mty::Insertion> insertions;

        std::vector<csl::Tensor>    momenta;

        std::vector<size_t> indices;

        std::vector<csl::Expr>      squaredMomenta;
    };

} // namespace mty