LorentzRSSpinor.h

// -*- C++ -*-
//
// LorentzRSSpinor.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_LorentzRSSpinor_H
#define ThePEG_LorentzRSSpinor_H

#include "ThePEG/Config/ThePEG.h"
#include "ThePEG/Vectors/ThreeVector.h"
#include "HelicityDefinitions.h"
#include "LorentzRSSpinor.fh"
#include "LorentzRSSpinorBar.h"
#include "LorentzSpinorBar.h"
#include "LorentzSpinor.h"
#include "LorentzPolarizationVector.h"

namespace ThePEG{
namespace Helicity{

template<typename Value>
class LorentzRSSpinor {

public:

  LorentzRSSpinor(SpinorType t = SpinorType::unknown) : _type(t), _spin() {}

  LorentzRSSpinor(complex<Value> a1, complex<Value> b1,
                  complex<Value> c1, complex<Value> d1,
                  complex<Value> a2, complex<Value> b2,
                  complex<Value> c2, complex<Value> d2,
                  complex<Value> a3, complex<Value> b3,
                  complex<Value> c3, complex<Value> d3,
                  complex<Value> a4, complex<Value> b4,
                  complex<Value> c4, complex<Value> d4,
                  SpinorType t=SpinorType::unknown) 
    : _type(t), _spin{{ {{a1,b1,c1,d1}},
                        {{a2,b2,c2,d2}},
                        {{a3,b3,c3,d3}},
                        {{a4,b4,c4,d4}}
                      }} {}
  
  template <typename U>
  LorentzRSSpinor(const LorentzRSSpinor<U> & other)
    : _type(other._type), _spin(other._spin) {}

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

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

  complex<Value> xs1() const {return _spin[0][0];}

  complex<Value> xs2() const {return _spin[0][1];}

  complex<Value> xs3() const {return _spin[0][2];}

  complex<Value> xs4() const {return _spin[0][3];}

  complex<Value> ys1() const {return _spin[1][0];}

  complex<Value> ys2() const {return _spin[1][1];}
  
  complex<Value> ys3() const {return _spin[1][2];}
  
  complex<Value> ys4() const {return _spin[1][3];}
  
  complex<Value> zs1() const {return _spin[2][0];}
  
  complex<Value> zs2() const {return _spin[2][1];}
  
  complex<Value> zs3() const {return _spin[2][2];}
  
  complex<Value> zs4() const {return _spin[2][3];}
  
  complex<Value> ts1() const {return _spin[3][0];}
  
  complex<Value> ts2() const {return _spin[3][1];}
  
  complex<Value> ts3() const {return _spin[3][2];}
  
  complex<Value> ts4() const {return _spin[3][3];}
  
  void setXS1(complex<Value> in) {_spin[0][0]=in;}
  
  void setXS2(complex<Value> in) {_spin[0][1]=in;}
  
  void setXS3(complex<Value> in) {_spin[0][2]=in;}
  
  void setXS4(complex<Value> in) {_spin[0][3]=in;}
  
  void setYS1(complex<Value> in) {_spin[1][0]=in;}
  
  void setYS2(complex<Value> in) {_spin[1][1]=in;}
  
  void setYS3(complex<Value> in) {_spin[1][2]=in;}
  
  void setYS4(complex<Value> in) {_spin[1][3]=in;}
  
  void setZS1(complex<Value> in) {_spin[2][0]=in;}
  
  void setZS2(complex<Value> in) {_spin[2][1]=in;}
  
  void setZS3(complex<Value> in) {_spin[2][2]=in;}
  
  void setZS4(complex<Value> in) {_spin[2][3]=in;}
  
  void setTS1(complex<Value> in) {_spin[3][0]=in;}
  
  void setTS2(complex<Value> in) {_spin[3][1]=in;}
  
  void setTS3(complex<Value> in) {_spin[3][2]=in;}
  
  void setTS4(complex<Value> in) {_spin[3][3]=in;}


  template <typename ValueB>
  LorentzRSSpinor<Value> & operator+=(const LorentzRSSpinor<ValueB> & a) {
    for(unsigned int ix=0;ix<4;++ix)
      for(unsigned int iy=0;iy<4;++iy)
        _spin[ix][iy] += a._spin[ix][iy];
    return *this;
  }
  
  template <typename ValueB>
  LorentzRSSpinor<Value> & operator-=(const LorentzRSSpinor<ValueB> & a) {
    for(unsigned int ix=0;ix<4;++ix)
      for(unsigned int iy=0;iy<4;++iy)
        _spin[ix][iy] -= a._spin[ix][iy];
    return *this;
  }
  
  LorentzRSSpinor<Value> & operator*=(double a) {
    for(unsigned int ix=0;ix<4;++ix)
      for(unsigned int iy=0;iy<4;++iy)
        _spin[ix][iy] *=a;
    return *this;
  }
  
  LorentzRSSpinor<Value> & operator/=(double a) {
    for(unsigned int ix=0;ix<4;++ix)
      for(unsigned int iy=0;iy<4;++iy)
        _spin[ix][iy] /=a;
    return *this;
  }
  
  LorentzSpinor<Value> dot(const LorentzPolarizationVector & vec) const {
    LorentzSpinor<Value> output(_type);
    complex<Value> temp;
    unsigned int ix;
    for(ix=0;ix<4;++ix) {
      temp  = _spin[3][ix]*vec.t();
      temp -= _spin[0][ix]*vec.x();
      temp -= _spin[1][ix]*vec.y();
      temp -= _spin[2][ix]*vec.z();
      output[ix]=temp;
    }
    return output;
  }

  LorentzSpinor<Value> dot(const LorentzMomentum & invec) const {
    LorentzSpinor<Value> output(_type);
    complex<Value> temp;
    LorentzVector<double> vec = UnitRemoval::InvE * invec;
    unsigned int ix;
    for(ix=0;ix<4;++ix) {
      temp  = - ( _spin[0][ix]*vec.x() + _spin[1][ix]*vec.y()+
                  _spin[2][ix]*vec.z() ) +  _spin[3][ix]*vec.t();
      output[ix]=temp;
    }
    return output;
  }

  LorentzRSSpinorBar<Value> bar() const;

  LorentzRSSpinor & boost(double,double,double);

  LorentzRSSpinor & boost(const Boost &);

  LorentzRSSpinor & transform(const LorentzRotation &);


  SpinorType Type() const {return _type;}

  template <typename ValueB>
  auto generalScalar(LorentzRSSpinorBar<ValueB>& fbar, 
                     Complex left, Complex right) 
  -> decltype(left*fbar(3,0)*this->ts1())
  {
    decltype(left*fbar(3,0)*ts1()) output; 
    unsigned int iz;
    output = 
      left*(fbar(3,0)*_spin[3][0]+fbar(3,1)*_spin[3][1])
      +right*(fbar(3,2)*_spin[3][2]+fbar(3,3)*_spin[3][3]);
    for(iz=0;iz<3;++iz) {
      output -=
        left*(fbar(iz,0)*_spin[iz][0]+fbar(iz,1)*_spin[iz][1])
        +right*(fbar(iz,2)*_spin[iz][2]+fbar(iz,3)*_spin[iz][3]);
    }
    return output;
  }
  
  template <typename ValueB>
  auto generalCurrent(LorentzSpinorBar<ValueB>& fbar, 
                      Complex left, Complex right) 
  -> LorentzVector<decltype(left*fbar.s1()*this->ts1())>
  {
    typedef decltype(left*fbar.s1()*ts1()) ResultT;
    ResultT output[4];
    for(size_t iz=0;iz<4;++iz)
      output[iz]= left*(fbar.s1()*_spin[iz][0]+fbar.s2()*_spin[iz][1])
        +right*(fbar.s3()*_spin[iz][2]+fbar.s4()*_spin[iz][3]);
    return LorentzVector<ResultT>(output[0],output[1],output[2],output[3]);
  }
  
private:

  SpinorType _type;

  std::array<std::array<complex<Value>,4>,4> _spin;
};


template <typename Value>
inline LorentzRSSpinor<double>
operator/(const LorentzRSSpinor<Value> & v, Value a) {
  return LorentzRSSpinor<double>(v.xs1()/a, v.xs2()/a, v.xs3()/a, v.xs4()/a,
                               v.ys1()/a, v.ys2()/a, v.ys3()/a, v.ys4()/a,
                               v.zs1()/a, v.zs2()/a, v.zs3()/a, v.zs4()/a,
                               v.ts1()/a, v.ts2()/a, v.ts3()/a, v.ts4()/a,
                               v.Type());
}

inline LorentzRSSpinor<double>
operator/(const LorentzRSSpinor<double> & v, Complex a) {
  return LorentzRSSpinor<double>(v.xs1()/a, v.xs2()/a, v.xs3()/a, v.xs4()/a,
                                 v.ys1()/a, v.ys2()/a, v.ys3()/a, v.ys4()/a,
                                 v.zs1()/a, v.zs2()/a, v.zs3()/a, v.zs4()/a,
                                 v.ts1()/a, v.ts2()/a, v.ts3()/a, v.ts4()/a,
                                 v.Type());
}

template <typename Value>
inline LorentzRSSpinor<Value> operator-(const LorentzRSSpinor<Value> & v) {
  return LorentzRSSpinor<Value>(-v.xs1(),-v.xs2(),-v.xs3(),-v.xs4(),
                                -v.ys1(),-v.ys2(),-v.ys3(),-v.ys4(),
                                -v.zs1(),-v.zs2(),-v.zs3(),-v.zs4(),
                                -v.ts1(),-v.ts2(),-v.ts3(),-v.ts4(),
                                v.Type());
}

template <typename ValueA, typename ValueB>
inline LorentzRSSpinor<ValueA>
operator+(LorentzRSSpinor<ValueA> a, const LorentzRSSpinor<ValueB> & b) {
  return a += b;
}

template <typename ValueA, typename ValueB>
inline LorentzRSSpinor<ValueA>
operator-(LorentzRSSpinor<ValueA> a, const LorentzRSSpinor<ValueB> & b) {
  return a -= b;
}

template <typename Value>
inline LorentzRSSpinor<Value>
operator*(const LorentzRSSpinor<Value> & a, double b) {
  return LorentzRSSpinor<Value>(a.xs1()*b, a.xs2()*b, a.xs3()*b, a.xs4()*b,
                                a.ys1()*b, a.ys2()*b, a.ys3()*b, a.ys4()*b,
                                a.zs1()*b, a.zs2()*b, a.zs3()*b, a.zs4()*b,
                                a.ts1()*b, a.ts2()*b, a.ts3()*b, a.ts4()*b,a.Type());
}

template <typename Value>
inline LorentzRSSpinor<Value>
operator*(double b, LorentzRSSpinor<Value> a) {
  return a *= b;
}
  
template <typename Value>
inline LorentzRSSpinor<Value>
operator*(const LorentzRSSpinor<Value> & a, Complex b) {
  return LorentzRSSpinor<Value>(a.xs1()*b, a.xs2()*b, a.xs3()*b, a.xs4()*b,
                                a.ys1()*b, a.ys2()*b, a.ys3()*b, a.ys4()*b,
                                a.zs1()*b, a.zs2()*b, a.zs3()*b, a.zs4()*b,
                                a.ts1()*b, a.ts2()*b, a.ts3()*b, a.ts4()*b,a.Type());
}

template <typename ValueA, typename ValueB>
inline auto operator*(complex<ValueB> a, const LorentzRSSpinor<ValueA> & v) 
  -> LorentzRSSpinor<decltype(a.real()*v.xs1().real())>
{
  return {a*v.xs1(), a*v.xs2(), a*v.xs3(), a*v.xs4(),
          a*v.ys1(), a*v.ys2(), a*v.ys3(), a*v.ys4(),
          a*v.zs1(), a*v.zs2(), a*v.zs3(), a*v.zs4(),
          a*v.ts1(), a*v.ts2(), a*v.ts3(), a*v.ts4(),v.Type()};
}

template <typename ValueA, typename ValueB>
inline auto operator*(const LorentzRSSpinor<ValueA> & v, complex<ValueB> b) 
  -> LorentzRSSpinor<decltype(b.real()*v.xs1().real())>
{
  return b*v;
}

template <typename ValueA, typename ValueB>
inline auto operator/(const LorentzRSSpinor<ValueA> & v, complex<ValueB> b) 
  -> LorentzRSSpinor<decltype(v.xs1().real()/b.real())>
{
  return {v.xs1()/b, v.xs2()/b, v.xs3()/b, v.xs4()/b,
          v.ys1()/b, v.ys2()/b, v.ys3()/b, v.ys4()/b,
          v.zs1()/b, v.zs2()/b, v.zs3()/b, v.zs4()/b,
          v.ts1()/b, v.ts2()/b, v.ts3()/b, v.ts4()/b,v.Type()};
}
  
template <typename ValueA, typename ValueB>
inline auto operator*(ValueB a, const LorentzRSSpinor<ValueA> & v) 
  -> LorentzRSSpinor<decltype(a*v.xs1().real())>
{
  return {a*v.xs1(), a*v.xs2(), a*v.xs3(), a*v.xs4(),
          a*v.ys1(), a*v.ys2(), a*v.ys3(), a*v.ys4(),
          a*v.zs1(), a*v.zs2(), a*v.zs3(), a*v.zs4(),
          a*v.ts1(), a*v.ts2(), a*v.ts3(), a*v.ts4(),v.Type()};
}

template <typename ValueA, typename ValueB>
inline auto operator*(const LorentzRSSpinor<ValueA> & v, ValueB b) 
  -> LorentzRSSpinor<decltype(b*v.xs1().real())>
{
  return b*v;
}

template <typename ValueA, typename ValueB>
inline auto operator/(const LorentzRSSpinor<ValueA> & v, ValueB b) 
  -> LorentzRSSpinor<decltype(v.xs1().real()/b)>
{
  return {v.xs1()/b, v.xs2()/b, v.xs3()/b, v.xs4()/b,
          v.ys1()/b, v.ys2()/b, v.ys3()/b, v.ys4()/b,
          v.zs1()/b, v.zs2()/b, v.zs3()/b, v.zs4()/b,
          v.ts1()/b, v.ts2()/b, v.ts3()/b, v.ts4()/b,v.Type()};
}

}
}
#ifndef ThePEG_TEMPLATES_IN_CC_FILE
#include "LorentzRSSpinor.tcc"
#endif 

#endif