LorentzMatrix.hh

Go to the documentation of this file.

#pragma once

#include "boca/math/LorentzVectorBase.hh"
#include "boca/math/Matrix3.hh"

namespace boca
{

template<typename Value_>
class LorentzMatrix : public boost::totally_ordered<LorentzMatrix<Value_>>
            , boost::additive<LorentzMatrix<Value_>>
{

    template<typename Value_2_>
    using ValueProduct = ValueProduct<Value_, Value_2_>;
    using ValueSquare = boca::ValueSquare<Value_>;
    using ValueCubed = boca::ValueCubed<Value_>;
    using Value4 = boca::Value4<Value_>;
    template<typename Value_2>
    using ValueQuotient = ValueQuotient<Value_, Value_2>;
    using ValueInverse = boost::units::divide_typeof_helper<double, Value_>;
    template<typename Value_2>
    using OnlyIfNotOrSameQuantity = typename std::enable_if < !IsQuantity<Value_2>::value || std::is_same<Value_, Value_2>::value >::type;
    template<typename Value_2>
    using OnlyIfNotQuantity = typename std::enable_if < !IsQuantity<Value_2>::value >::type;

public:

    LorentzMatrix() {}

    LorentzMatrix(Matrix3<Value_> const &matrix, Vector3<Value_> const &vector, Value_ scalar, Matrix symmetry = Matrix::none)
    {
        if (symmetry == Matrix::none) symmetry == matrix.Symmetry();
        if (symmetry != Matrix::symmetric || symmetry != Matrix::antisymmetric) std::cout << "matrix symmetry not well defiend" << '\n';
        SetMatrix(matrix);
        SetVector(vector, symmetry);
        SetScalar(scalar);
    }

    LorentzMatrix(Matrix3<Value_> const &matrix)
    {
        SetMatrix(matrix);
        SetVector({0, 0, 0}, Matrix::symmetric);
        SetScalar(Value_(1));
    }

    LorentzMatrix(Vector3<double> const &vector)
    {
        SetBoost(vector);
    }

    LorentzMatrix(LorentzVectorBase<Value_> const &x, LorentzVectorBase<Value_> const &y, LorentzVectorBase<Value_> const &z, LorentzVectorBase<Value_> const &t, Matrix matrix = Matrix::row)
    {
        switch (matrix) {
        case Matrix::row:
            SetRows(x, y, z, t);
            break;
        case Matrix::column :
            SetColumns(x, y, z, t);
            break;
        default :
            std::cout << "Maformed matrix constructor: " << Name(matrix) << '\n';
        }
    }


    void SetMatrix(Matrix3<Value_> const &matrix)
    {
        x_.SetVector(matrix.X());
        y_.SetVector(matrix.Y());
        z_.SetVector(matrix.Z());
    }

    void SetVector(Vector3<Value_> const &vector, Matrix symmetry)
    {
        x_.SetScalar(vector.X());
        y_.SetScalar(vector.Y());
        z_.SetScalar(vector.Z());
        switch (symmetry) {
        case Matrix::symmetric :
            t_.SetVector(vector);
            break;
        case Matrix::antisymmetric :
            t_.SetVector(-vector);
            break;
        default :
            std::cout << "matrix symmetry not well defiend" << '\n';
        }
    }

    void SetScalar(Value_ scalar)
    {
        t_.SetScalar(scalar);
    }

    void SetRows(LorentzVectorBase<Value_> const &x, LorentzVectorBase<Value_> const &y, LorentzVectorBase<Value_> const &z, LorentzVectorBase<Value_> const &t)
    {
        x_ = x;
        y_ = y;
        z_ = z;
        t_ = t;
    }

    void SetColumns(LorentzVectorBase<Value_> const &x, LorentzVectorBase<Value_> const &y, LorentzVectorBase<Value_> const &z, LorentzVectorBase<Value_> const &t)
    {
        x_ = {x.X(), y.X(), z.X(), t.X()};
        y_ = {x.Y(), y.Y(), z.Y(), t.Y()};
        z_ = {x.Z(), y.Z(), z.Z(), t.Z()};
        t_ = {x.T(), y.T(), z.T(), t.T()};
    }

    LorentzVectorBase<Value_> const& X() const
    {
        return x_;
    }

    LorentzVectorBase<Value_> &X()
    {
        return x_;
    }

    LorentzVectorBase<Value_> const& Y() const
    {
        return y_;
    }

    LorentzVectorBase<Value_> &Y()
    {
        return y_;
    }

    LorentzVectorBase<Value_> const& Z() const
    {
        return z_;
    }

    LorentzVectorBase<Value_> &Z()
    {
        return z_;
    }

    LorentzVectorBase<Value_> T() const
    {
        return t_;
    }

    LorentzVectorBase<Value_> &T()
    {
        return t_;
    }


    Value4 Determinant()const
    {
        return  boost::accumulate(LorentzDimensions(), Value4(0), [&](Value4 & sum, LorentzDim dim) {
            return sum + Laplace(LorentzDim::x, dim);
        });
    }

    Value4 Laplace(LorentzDim dim_1, LorentzDim dim_2) const
    {
        return (*this)[dim_1][dim_2] * Cofactor(dim_1, dim_2);
    }

    ValueCubed Cofactor(LorentzDim dim_1, LorentzDim dim_2) const
    {
        return static_cast<double>(Sign(dim_1, dim_2)) * Minor(dim_1, dim_2);
    }

    ValueCubed Minor(LorentzDim delete_1, LorentzDim delete_2) const
    {
        return SubMatrix(delete_1, delete_2).Determinant();
    }

    Value4 ReducedDeterminant(LorentzDim dim_1, LorentzDim dim_2) const
    {
        return Determinant() - (*this)[dim_1][dim_2] * Cofactor(dim_1, dim_2);
    }

    Value4 ReducedLaplace(LorentzDim dim_1, LorentzDim dim_2, Dim3 dim_3, Dim3 dim_4)
    {
        return (*this)[dim_1][dim_2] * SubMatrix(dim_1, dim_2).ReducedDeterminant(dim_3, dim_4);
    }

    int Sign(LorentzDim i, LorentzDim j) const
    {
        return (static_cast<int>(i) + static_cast<int>(j)) % 2 ? -1 : 1;
    }


    Matrix3<Value_> SubMatrix(LorentzDim delete_1, LorentzDim delete_2) const
    {
        auto dim3_1 = EnumIterator<Dim3> {Dim3::x};
        auto dim3_2 = EnumIterator<Dim3> {Dim3::x};
        auto matrix = Matrix3<Value_> {};
        for (auto dim4_1 : LorentzDimensions()) {
            if (dim4_1 == delete_1) continue;
            for (auto dim4_2 : LorentzDimensions()) {
                if (dim4_2 == delete_2) continue;
                matrix(*dim3_1, *dim3_2) = (*this)(dim4_1, dim4_2);
                ++dim3_2;
            }
            ++dim3_1;
            dim3_2.Set(Dim3::x);
        }
        return matrix;
    }


    LorentzVectorBase<Value_> ColumnX() const
    {
        return {x_.X(), y_.X(), z_.X(), t_.X()};
    }

    LorentzVectorBase<Value_> ColumnY() const
    {
        return {x_.Y(), y_.Y(), z_.Y(), t_.Y()};
    }

    LorentzVectorBase<Value_> ColumnZ() const
    {
        return {x_.Z(), y_.Z(), z_.Z(), t_.Z()};
    }

    LorentzVectorBase<Value_> ColumnT() const
    {
        return {x_.T(), y_.T(), z_.T(), t_.T()};
    }


    template<typename Value_2_>
    auto Multiply(const LorentzVectorBase<Value_2_> &vector) const
    {
        return {x_ * vector, y_ * vector, z_ * vector, t_ * vector};
    }

    template<typename Value_2_>
    LorentzMatrix Multiply(const LorentzMatrix<Value_2_> &matrix) const
    {
        return {{x_ * matrix.ColumnX(), x_ * matrix.ColumnY(), x_ * matrix.ColumnZ(), x_ * matrix.ColumnT()},
            {y_ * matrix.ColumnX(), y_ * matrix.ColumnY(), y_ * matrix.ColumnZ(), y_ * matrix.ColumnT()},
            {z_ * matrix.ColumnX(), z_ * matrix.ColumnY(), z_ * matrix.ColumnZ(), z_ * matrix.ColumnT()},
            {t_ * matrix.ColumnX(), t_ * matrix.ColumnY(), t_ * matrix.ColumnZ(), t_ * matrix.ColumnT()}
        };
    }

    template<typename Value_2_>
    LorentzMatrix &Transform(const LorentzMatrix<Value_2_> &matrix)
    {
        return *this = matrix.Multiply(*this);
    }

    template<typename Value_2_>
    auto &Transform(const Matrix3<Value_2_> &matrix)
    {
        return Transform(LorentzMatrix(matrix));
    }

    auto Transposed() const
    {
        return {ColumnX(), ColumnY(), ColumnZ(),ColumnT()};
    }

    auto &Transpose()
    {
        return *this = Transposed();
    }

    template<typename Value_2_>
    auto &Boost(const Vector3<Value_> &boost)
    {
        return Transform(LorentzMatrix(boost));
    }

    template<typename Value_2_>
    auto &Rotate(Angle angle, Dim3 dim)
    {
        return Transform(Matrix3<double>().Rotate(angle,  dim));
    }

    template<typename Value_2_>
    auto &Rotate(Angle angle, const Vector3<Value_> &axis)
    {
        return Transform(Matrix3<double>().Rotate(angle, axis));
    }


    bool operator<(LorentzMatrix const &matrix) const
    {
        return abs(Determinant()) < abs(matrix.Determinant());
    }

    bool operator==(LorentzMatrix const &matrix) const
    {
        return x_ == matrix.x_ && y_ == matrix.y_ && z_ == matrix.z_ && t_ == matrix.t_ ;
    }

    template<typename Value_2_>
    auto operator*(const LorentzVectorBase<Value_2_> &vector) const
    {
        return Multiply(vector);
    }

    template<typename Value_2_>
    auto operator*(const LorentzMatrix<Value_2_> &matrix) const
    {
        return Multiply(matrix);
    }

    auto &operator*= (const LorentzMatrix<Value_2_> &matrix)
    {
        return *this = Multiply(matrix);
    }

    LorentzVectorBase<Value_> const &operator[](LorentzDim lorentz_dim) const
    {
        switch (lorentz_dim) {
        case LorentzDim::x :
            return x_;
        case LorentzDim::y :
            return y_;
        case LorentzDim::z :
            return z_;
        case LorentzDim::t :
            return t_;
        default : Default("LorentzMatrix",  Name(lorentz_dim));
        return x_;
        }
    }

    LorentzVectorBase<Value_> &operator[](LorentzDim lorentz_dim)
    {
        return const_cast<LorentzVectorBase<Value_> &>(static_cast<LorentzMatrix<Value_> const &>(*this)[lorentz_dim]);
    }

    friend auto &operator<<(std::ostream &stream, LorentzMatrix<Value_> const &matrix)
    {
        for(auto const& vector : matrix) stream << vector;
        return stream;
    }


private:

    LorentzVectorBase<Value_> x_;

    LorentzVectorBase<Value_> y_;

    LorentzVectorBase<Value_> z_;

    LorentzVectorBase<Value_> t_;

    // Set elements according to a boost vector.
    //boost this Lorentz vector
    void SetBoost(Vector3<double> const &boost)
    {
        auto bp2 = boost.Mag2();
        auto gamma = 1 / std::sqrt(1 - bp2);
        auto bgamma = sqr(gamma) / (1 + gamma);
        auto matrix = Matrix3<double>(bgamma, Matrix::uniform) + Multiply(boost, boost) + Matrix3<double>(1);
        auto vector = boost + Vector3<double>(gamma, gamma, gamma);
        *this = LorentzMatrix<double>(matrix, vector, gamma);
    }


    // Assignment.
//     LorentzMatrix<Value_> &operator=(const Matrix3<Value_> &matrix)
//     {
//         SetMatrix(matrix);
//         SetVector({0, 0, 0}, Matrix::symmetric);
//         SetScalar(Value_(1));
//         return *this;
//     }


};

}