LorentzTensor.h

// -*- C++ -*-
//
// LorentzTensor.h is a part of ThePEG - Toolkit for HEP Event Generation
// Copyright (C) 2003-2019 Peter Richardson, Leif Lonnblad
//
// ThePEG is licenced under version 3 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
#ifndef ThePEG_LorentzTensor_H
#define ThePEG_LorentzTensor_H
// This is the declaration of the LorentzTensor class.

#include "ThePEG/Config/PhysicalQtyComplex.h"
#include "ThePEG/Config/ThePEG.h"
#include "LorentzPolarizationVector.h"

namespace ThePEG {
namespace Helicity {

// compiler magic needs these pre-declarations to make friend templates work
template<typename Value> class LorentzTensor;

template<typename Value> 
class LorentzTensor {

public:

  LorentzTensor() = default;

  LorentzTensor(complex<Value> xx, complex<Value> xy, 
                complex<Value> xz, complex<Value> xt,
                complex<Value> yx, complex<Value> yy,
                complex<Value> yz, complex<Value> yt,
                complex<Value> zx, complex<Value> zy,
                complex<Value> zz, complex<Value> zt,
                complex<Value> tx, complex<Value> ty,
                complex<Value> tz, complex<Value> tt)
  : _tensor{{ {{xx,xy,xz,xt}},
              {{yx,yy,yz,yt}},
              {{zx,zy,zz,zt}},
              {{tx,ty,tz,tt}} }} {}
  LorentzTensor(const LorentzPolarizationVector & p,
                const LorentzPolarizationVector & q) {
    setXX(p.x() * q.x()); setYX(p.y() * q.x());
    setZX(p.z() * q.x()); setTX(p.t() * q.x());
    setXY(p.x() * q.y()); setYY(p.y() * q.y());
    setZY(p.z() * q.y()); setTY(p.t() * q.y());
    setXZ(p.x() * q.z()); setYZ(p.y() * q.z());
    setZZ(p.z() * q.z()); setTZ(p.t() * q.z());
    setXT(p.x() * q.t()); setYT(p.y() * q.t());
    setZT(p.z() * q.t()); setTT(p.t() * q.t());
  }

  complex<Value> xx() const {return _tensor[0][0];}

  complex<Value> yx() const {return _tensor[1][0];}
  complex<Value> zx() const {return _tensor[2][0];}

  complex<Value> tx() const {return _tensor[3][0];}

  complex<Value> xy() const {return _tensor[0][1];}

  complex<Value> yy() const {return _tensor[1][1];}

  complex<Value> zy() const {return _tensor[2][1];}

  complex<Value> ty() const {return _tensor[3][1];}

  complex<Value> xz() const {return _tensor[0][2];}

  complex<Value> yz() const {return _tensor[1][2];}

  complex<Value> zz() const {return _tensor[2][2];}

  complex<Value> tz() const {return _tensor[3][2];}

  complex<Value> xt() const {return _tensor[0][3];}

  complex<Value> yt() const {return _tensor[1][3];}

  complex<Value> zt() const {return _tensor[2][3];}

  complex<Value> tt() const {return _tensor[3][3];}

  void setXX(complex<Value> a) {_tensor[0][0]=a;}

  void setYX(complex<Value> a) {_tensor[1][0]=a;}

  void setZX(complex<Value> a) {_tensor[2][0]=a;}

  void setTX(complex<Value> a) {_tensor[3][0]=a;}

  void setXY(complex<Value> a) {_tensor[0][1]=a;}

  void setYY(complex<Value> a) {_tensor[1][1]=a;}

  void setZY(complex<Value> a) {_tensor[2][1]=a;}

  void setTY(complex<Value> a) {_tensor[3][1]=a;}

  void setXZ(complex<Value> a) {_tensor[0][2]=a;}

  void setYZ(complex<Value> a) {_tensor[1][2]=a;}

  void setZZ(complex<Value> a) {_tensor[2][2]=a;}

  void setTZ(complex<Value> a) {_tensor[3][2]=a;}

  void setXT(complex<Value> a) {_tensor[0][3]=a;}

  void setYT(complex<Value> a) {_tensor[1][3]=a;}

  void setZT(complex<Value> a) {_tensor[2][3]=a;}

  void setTT(complex<Value> a) {_tensor[3][3]=a;}

  complex<Value> operator () (int i, int j) const {
    assert( i>=0 && i<=3 && j>=0 && j<=3);
    return _tensor[i][j];
  }

  complex<Value> & operator () (int i, int j) {
    assert( i>=0 && i<=3 && j>=0 && j<=3);
    return _tensor[i][j];
  }

  LorentzTensor & boost(double,double,double);

  LorentzTensor<Value> & boost(const Boost & b) {
    return boost(b.x(), b.y(), b.z());
  }

  LorentzTensor & transform(const SpinOneLorentzRotation & r){
    unsigned int ix,iy,ixa,iya;
    LorentzTensor<Value> output;
    complex<Value> temp;
    for(ix=0;ix<4;++ix) {
      for(iy=0;iy<4;++iy) {
        temp=complex<Value>();
        for(ixa=0;ixa<4;++ixa) {
          for(iya=0;iya<4;++iya)
            temp+=r(ix,ixa)*r(iy,iya)*(*this)(ixa,iya);
        }
        output(ix,iy)=temp;
      }
    }
    *this=output;
    return *this;
  }
  
  LorentzTensor<Value> conjugate() {
    return LorentzTensor<Value>(conj(xx()), conj(xy()), conj(xz()), conj(xt()),
                                conj(yx()), conj(yy()), conj(yz()), conj(yt()),
                                conj(zx()), conj(zy()), conj(zz()), conj(zt()),
                                conj(tx()), conj(ty()), conj(tz()), conj(tt()));
  }


  LorentzTensor<Value> operator*=(Complex a) {
    for(int ix=0;ix<4;++ix)
      for(int iy=0;iy<4;++iy) _tensor[ix][iy]*=a;
    return *this;
  }

  template <typename T, typename U>
  friend auto
  operator*(const LorentzTensor<T> & t, const LorentzTensor<U> & u)
  -> decltype(t.xx()*u.xx())
  ;
    
  LorentzTensor<Value> operator+(const LorentzTensor<Value> & in) const {
    return LorentzTensor<Value>(xx()+in.xx(),xy()+in.xy(),xz()+in.xz(),xt()+in.xt(),
                                yx()+in.yx(),yy()+in.yy(),yz()+in.yz(),yt()+in.yt(),
                                zx()+in.zx(),zy()+in.zy(),zz()+in.zz(),zt()+in.zt(),
                                tx()+in.tx(),ty()+in.ty(),tz()+in.tz(),tt()+in.tt());
  }
  
  LorentzTensor<Value> operator-(const LorentzTensor<Value> & in) const {
    return LorentzTensor<Value>(xx()-in.xx(),xy()-in.xy(),xz()-in.xz(),xt()-in.xt(),
                                yx()-in.yx(),yy()-in.yy(),yz()-in.yz(),yt()-in.yt(),
                                zx()-in.zx(),zy()-in.zy(),zz()-in.zz(),zt()-in.zt(),
                                tx()-in.tx(),ty()-in.ty(),tz()-in.tz(),tt()-in.tt());
  }

  complex<Value> trace() const {
    return _tensor[3][3]-_tensor[0][0]-_tensor[1][1]-_tensor[2][2];
  }

  LorentzVector<complex<Value> > preDot (const LorentzPolarizationVector & vec) const {
    LorentzVector<complex<Value> > output;
    output.setX(vec.t()*_tensor[3][0]-vec.x()*_tensor[0][0]-
                vec.y()*_tensor[1][0]-vec.z()*_tensor[2][0]);
    output.setY(vec.t()*_tensor[3][1]-vec.x()*_tensor[0][1]-
                vec.y()*_tensor[1][1]-vec.z()*_tensor[2][1]);
    output.setZ(vec.t()*_tensor[3][2]-vec.x()*_tensor[0][2]-
                vec.y()*_tensor[1][2]-vec.z()*_tensor[2][2]);
    output.setT(vec.t()*_tensor[3][3]-vec.x()*_tensor[0][3]-
                vec.y()*_tensor[1][3]-vec.z()*_tensor[2][3]);
    return output;
  }

  LorentzVector<complex<Value> > postDot(const LorentzPolarizationVector & vec) const {
    LorentzVector<complex<Value> > output;
    output.setX(vec.t()*_tensor[0][3]-vec.x()*_tensor[0][0]-
                vec.y()*_tensor[0][1]-vec.z()*_tensor[0][2]);
    output.setY(vec.t()*_tensor[1][3]-vec.x()*_tensor[1][0]-
                vec.y()*_tensor[1][1]-vec.z()*_tensor[1][2]);
    output.setZ(vec.t()*_tensor[2][3]-vec.x()*_tensor[2][0]-
                vec.y()*_tensor[2][1]-vec.z()*_tensor[2][2]);
    output.setT(vec.t()*_tensor[3][3]-vec.x()*_tensor[3][0]-
                vec.y()*_tensor[3][1]-vec.z()*_tensor[3][2]);
    return output;
  }

  auto preDot (const Lorentz5Momentum & vec) const 
  -> LorentzVector<decltype(vec.x()*this->xx())>
  {
    LorentzVector<decltype(vec.x()*this->xx())> output;
    output.setX(vec.t()*_tensor[3][0]-vec.x()*_tensor[0][0]-
                vec.y()*_tensor[1][0]-vec.z()*_tensor[2][0]);
    output.setY(vec.t()*_tensor[3][1]-vec.x()*_tensor[0][1]-
                vec.y()*_tensor[1][1]-vec.z()*_tensor[2][1]);
    output.setZ(vec.t()*_tensor[3][2]-vec.x()*_tensor[0][2]-
                vec.y()*_tensor[1][2]-vec.z()*_tensor[2][2]);
    output.setT(vec.t()*_tensor[3][3]-vec.x()*_tensor[0][3]-
                vec.y()*_tensor[1][3]-vec.z()*_tensor[2][3]);
    return output;
  }

  auto postDot(const Lorentz5Momentum & vec) const 
  -> LorentzVector<decltype(vec.x()*this->xx())>
  {
    LorentzVector<decltype(vec.x()*this->xx())> output;
    output.setX(vec.t()*_tensor[0][3]-vec.x()*_tensor[0][0]-
                vec.y()*_tensor[0][1]-vec.z()*_tensor[0][2]);
    output.setY(vec.t()*_tensor[1][3]-vec.x()*_tensor[1][0]-
                vec.y()*_tensor[1][1]-vec.z()*_tensor[1][2]);
    output.setZ(vec.t()*_tensor[2][3]-vec.x()*_tensor[2][0]-
                vec.y()*_tensor[2][1]-vec.z()*_tensor[2][2]);
    output.setT(vec.t()*_tensor[3][3]-vec.x()*_tensor[3][0]-
                vec.y()*_tensor[3][1]-vec.z()*_tensor[3][2]);
    return output;
  }
private:

  std::array<std::array<complex<Value>,4>,4> _tensor;

};

template<typename T, typename U> 
inline auto
operator*(complex<U> a, const LorentzTensor<T> & t) 
-> LorentzTensor<decltype(a.real()*t.xx().real())>
{
  return 
    {a*t.xx(), a*t.xy(), a*t.xz(), a*t.xt(),
     a*t.yx(), a*t.yy(), a*t.yz(), a*t.yt(),
     a*t.zx(), a*t.zy(), a*t.zz(), a*t.zt(),
     a*t.tx(), a*t.ty(), a*t.tz(), a*t.tt()};
}

template<typename T, typename U> 
inline auto
operator*(const LorentzTensor<T> & t,complex<U> a) 
-> LorentzTensor<decltype(a.real()*t.xx().real())>
{
  return 
    {a*t.xx(), a*t.xy(), a*t.xz(), a*t.xt(),
     a*t.yx(), a*t.yy(), a*t.yz(), a*t.yt(),
     a*t.zx(), a*t.zy(), a*t.zz(), a*t.zt(),
     a*t.tx(), a*t.ty(), a*t.tz(), a*t.tt()};
}

template<typename T, typename U> 
inline auto
operator*(const LorentzVector<U> & v, 
          const LorentzTensor<T> & t) 
-> LorentzVector<decltype(v.t()*t(3,0))>
{
  LorentzVector<decltype(v.t()*t(3,0))> outvec;
  outvec.setX( v.t()*t(3,0)-v.x()*t(0,0)
              -v.y()*t(1,0)-v.z()*t(2,0));
  outvec.setY( v.t()*t(3,1)-v.x()*t(0,1)
              -v.y()*t(1,1)-v.z()*t(2,1));
  outvec.setZ( v.t()*t(3,2)-v.x()*t(0,2)
              -v.y()*t(1,2)-v.z()*t(2,2));
  outvec.setT( v.t()*t(3,3)-v.x()*t(0,3)
              -v.y()*t(1,3)-v.z()*t(2,3));
  return outvec;
}

template<typename T, typename U> 
inline auto
operator*(const LorentzTensor<T> & t, const LorentzVector<U> & v)
-> LorentzVector<decltype(v.t()*t(0,3))>
{
  LorentzVector<decltype(v.t()*t(0,3))> outvec;
  outvec.setX( v.t()*t(0,3)-v.x()*t(0,0)
              -v.y()*t(0,1)-v.z()*t(0,2));
  outvec.setY( v.t()*t(1,3)-v.x()*t(1,0)
              -v.y()*t(1,1)-v.z()*t(1,2));
  outvec.setZ( v.t()*t(2,3)-v.x()*t(2,0)
              -v.y()*t(2,1)-v.z()*t(2,2));
  outvec.setT( v.t()*t(3,3)-v.x()*t(3,0)
              -v.y()*t(3,1)-v.z()*t(3,2));
  return outvec;
}

template <typename T, typename U>
inline auto
operator*(const LorentzTensor<T> & t, 
          const LorentzTensor<U> & u) 
-> decltype(t.xx()*u.xx())
{
  using RetT = decltype(t.xx()*u.xx());
  RetT output=RetT(),temp;
  for(unsigned int ix=0;ix<4;++ix) {
    temp = t._tensor[ix][3]*u._tensor[ix][3];
    for(unsigned int iy=0;iy<3;++iy) {
      temp+= t._tensor[ix][iy]*u._tensor[ix][iy];
    }
    if(ix<3) output-=temp;
    else     output+=temp;
  }
  return output;
}

}
}

#ifndef ThePEG_TEMPLATES_IN_CC_FILE
#include "LorentzTensor.tcc"
#endif 

#endif