LorentzSpinor.h

// -*- C++ -*-
//
// LorentzSpinor.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_LorentzSpinor_H
#define ThePEG_LorentzSpinor_H
// This is the declaration of the LorentzSpinor class.
#include "ThePEG/Config/ThePEG.h"
#include "ThePEG/Vectors/LorentzRotation.h"
#include "ThePEG/Vectors/ThreeVector.h"
#include "HelicityDefinitions.h"
#include "LorentzSpinor.fh"
#include "LorentzSpinorBar.h"
#include "LorentzPolarizationVector.h"
#include "LorentzTensor.h"
#include <array>
 
namespace ThePEG{
namespace Helicity{
 
template<typename Value>
class LorentzSpinor {
public:
 
  LorentzSpinor(SpinorType t = SpinorType::unknown) : _type(t), _spin() {}
 
  LorentzSpinor(complex<Value> a,complex<Value> b,
                complex<Value> c,complex<Value> d,
                SpinorType s = SpinorType::unknown) : _type(s), _spin{{a,b,c,d}} {}
 
  template <typename U>
  LorentzSpinor(const LorentzSpinor<U> & other)
    : _type(other._type), _spin(other._spin) {}
 
  complex<Value> operator[](int i) const  {
    assert( i >= 0 && i <= 3 );
    return _spin[i];
  }
 
  complex<Value> operator()(int i) const  {
    assert( i >= 0 && i <= 3 );
    return _spin[i];
  }
  
  complex<Value> & operator()(int i) {
    assert( i >= 0 && i <= 3 );
    return _spin[i];
  }
  
  complex<Value> & operator[](int i) {
    assert( i >= 0 && i <= 3 );
    return _spin[i];
  }
  
  complex<Value> s1() const {return _spin[0];}
 
  complex<Value> s2() const {return _spin[1];}
 
  complex<Value> s3() const {return _spin[2];}
 
  complex<Value> s4() const {return _spin[3];}
 
  void setS1(complex<Value> in) {_spin[0]=in;}
 
  void setS2(complex<Value> in) {_spin[1]=in;}
 
  void setS3(complex<Value> in) {_spin[2]=in;}
 
  void setS4(complex<Value> in) {_spin[3]=in;}
 
 
  template <typename ValueB>
  LorentzSpinor<Value> & operator+=(const LorentzSpinor<ValueB> & a) {
    for(unsigned int ix=0;ix<4;++ix) _spin[ix] += a._spin[ix];
    return *this;
  }
  
  template <typename ValueB>
  LorentzSpinor<Value> & operator-=(const LorentzSpinor<ValueB> & a) {
    for(unsigned int ix=0;ix<4;++ix) _spin[ix] -= a._spin[ix];
    return *this;
  }
  
  LorentzSpinor<Value> & operator*=(double a) {
    for(unsigned int ix=0;ix<4;++ix) _spin[ix] *=a;
    return *this;
  }
 
  LorentzSpinor<Value> & operator/=(double a) {
    for(unsigned int ix=0;ix<4;++ix) _spin[ix] /=a;
    return *this;
  }
  
  LorentzSpinorBar<Value> bar() const;
 
  LorentzSpinor conjugate() const;
 
  LorentzSpinor & boost(double,double,double);
 
  LorentzSpinor & boost(const Boost &);
 
  LorentzSpinor & transform(const SpinHalfLorentzRotation & );
 
  LorentzSpinor & transform(const LorentzRotation & r) {
    transform(r.half());
    return *this;
  }
 
  SpinorType Type() const {return _type;}
 
  template<typename ValueB> 
  auto projectionOperator(const LorentzVector<ValueB> & p, 
                          const ValueB & m) const 
  -> LorentzSpinor<decltype(m*Value())>
  {
    LorentzSpinor<decltype(m*Value())> spin;
    static const Complex ii(0.,1.);
    complex<ValueB> p0pp3=p.t()+p.z();
    complex<ValueB> p0mp3=p.t()-p.z();
    complex<ValueB> p1pp2=p.x()+ii*p.y();
    complex<ValueB> p1mp2=p.x()-ii*p.y();
    spin.setS1(m*s1()+p0mp3*s3()-p1mp2*s4());
    spin.setS2(m*s2()+p0pp3*s4()-p1pp2*s3());
    spin.setS3(m*s3()+p0pp3*s1()+p1mp2*s2());
    spin.setS4(m*s4()+p0mp3*s2()+p1pp2*s1());
    return spin;
  }
 
  LorentzSpinor
  helicityProjectionOperator(const Complex & gL, const Complex & gR) const {
    LorentzSpinor spin;
    spin.setS1(gL*s1());
    spin.setS2(gL*s2());
    spin.setS3(gR*s3());
    spin.setS4(gR*s4());
    return spin;
  }
 
 
  template<typename ValueB> 
  auto slash(const LorentzVector<ValueB> & p) const 
  -> LorentzSpinor<decltype(p.t()*Value())>
  {
    LorentzSpinor<decltype(p.t()*Value())> spin;
    static const Complex ii(0.,1.);
    complex<ValueB> p0pp3=p.t()+p.z();
    complex<ValueB> p0mp3=p.t()-p.z();
    complex<ValueB> p1pp2=p.x()+ii*p.y();
    complex<ValueB> p1mp2=p.x()-ii*p.y();
    spin.setS1(p0mp3*s3()-p1mp2*s4());
    spin.setS2(p0pp3*s4()-p1pp2*s3());
    spin.setS3(p0pp3*s1()+p1mp2*s2());
    spin.setS4(p0mp3*s2()+p1pp2*s1());
    return spin;
  }
 
  template<typename ValueB> 
  auto slash(const LorentzVector<complex<ValueB> > & p) const 
  -> LorentzSpinor<decltype(ValueB()*Value())>
  {
    LorentzSpinor<decltype(ValueB()*Value())> spin;
    static const Complex ii(0.,1.);
    complex<ValueB> p0pp3=p.t()+p.z();
    complex<ValueB> p0mp3=p.t()-p.z();
    complex<ValueB> p1pp2=p.x()+ii*p.y();
    complex<ValueB> p1mp2=p.x()-ii*p.y();
    spin.setS1(p0mp3*s3()-p1mp2*s4());
    spin.setS2(p0pp3*s4()-p1pp2*s3());
    spin.setS3(p0pp3*s1()+p1mp2*s2());
    spin.setS4(p0mp3*s2()+p1pp2*s1());
    return spin;
  }
 
  template<typename ValueB>
  auto leftCurrent(const LorentzSpinorBar<ValueB>& fb) const 
  -> LorentzVector<decltype(fb.s3()*this->s2())>
  {
    typedef decltype(fb.s3()*s2()) ResultT;
    LorentzVector<ResultT> vec;
    Complex ii(0.,1.);
    ResultT p1(fb.s3()*s2()),p2(fb.s4()*s1());
    vec.setX(   -(p1+p2) );
    vec.setY( ii*(p1-p2) );
    p1 = fb.s3()*s1();p2 = fb.s4()*s2();
    vec.setZ(   -(p1-p2) );
    vec.setT(    (p1+p2) );
    return vec;
  }
 
  template<typename ValueB>
  auto rightCurrent(const LorentzSpinorBar<ValueB>& fb) const 
  -> LorentzVector<decltype(fb.s1()*this->s4())>
  {
    typedef decltype(fb.s1()*s4()) ResultT;
    LorentzVector<ResultT> vec;
    Complex ii(0.,1.);
    ResultT p1(fb.s1()*s4()),p2(fb.s2()*s3());
    vec.setX(     (p1+p2));
    vec.setY( -ii*(p1-p2));
    p1 = fb.s1()*s3();p2 = fb.s2()*s4();
    vec.setZ(     (p1-p2));
    vec.setT(     (p1+p2));
    return vec;
  }
 
  template<typename ValueB>
  auto vectorCurrent(const LorentzSpinorBar<ValueB>& fb) const 
  -> LorentzVector<decltype(fb.s1()*this->s4())>
  {
    typedef decltype(fb.s1()*this->s4()) ResultT;
    LorentzVector<ResultT> vec;
    Complex ii(0.,1.);
    ResultT s1s4(fb.s1()*s4()),s2s3(fb.s2()*s3()),
      s3s2(fb.s3()*s2()),s4s1(fb.s4()*s1()),
      s1s3(fb.s1()*s3()),s2s4(fb.s2()*s4()),
      s3s1(fb.s3()*s1()),s4s2(fb.s4()*s2());
    vec.setX(      s1s4+s2s3-s3s2-s4s1 );
    vec.setY( -ii*(s1s4-s2s3-s3s2+s4s1));
    vec.setZ(      s1s3-s2s4-s3s1+s4s2 );
    vec.setT(      s1s3+s2s4+s3s1+s4s2);
    return vec;
  }
 
  template<typename ValueB>
  auto generalCurrent(const LorentzSpinorBar<ValueB>& fb,
                      Complex left, Complex right) const
  -> LorentzVector<decltype(fb.s3()*this->s2())>
  {
    typedef decltype(fb.s3()*this->s2()) ResultT;
    LorentzVector<ResultT> vec;
    Complex ii(0.,1.);
    ResultT p1(fb.s3()*s2()),p2(fb.s4()*s1());
    vec.setX(   -left*(p1+p2));
    vec.setY( ii*left*(p1-p2));
    p1 = fb.s3()*s1();p2 = fb.s4()*s2();
    vec.setZ(   -left*(p1-p2));
    vec.setT(    left*(p1+p2));
    p1=fb.s1()*s4();p2=fb.s2()*s3();
    vec.setX(vec.x()+right*(p1+p2));
    vec.setY(vec.y()-ii*right*(p1-p2));
    p1 = fb.s1()*s3();p2 = fb.s2()*s4();
    vec.setZ(vec.z()+right*(p1-p2));
    vec.setT(vec.t()+right*(p1+p2));
    return vec;
  }
 
  template<typename ValueB>
  auto leftScalar(const LorentzSpinorBar<ValueB>& fb) const  
  -> decltype(fb.s1()*this->s1())
  {
    return fb.s1()*s1()+fb.s2()*s2();
  }
 
  template<typename ValueB>
  auto rightScalar(const LorentzSpinorBar<ValueB>& fb) const 
  -> decltype(fb.s3()*this->s3())
  {
    return fb.s3()*s3()+fb.s4()*s4();
  }
  
  template<typename ValueB>
  auto scalar(const LorentzSpinorBar<ValueB>& fb) const 
  -> decltype(fb.s1()*this->s1())
  {
    return fb.s1()*s1()+fb.s2()*s2()
          +fb.s3()*s3()+fb.s4()*s4();
  }
 
  template<typename ValueB>
  auto pseudoScalar(const LorentzSpinorBar<ValueB>& fb) const 
  -> decltype(fb.s1()*this->s1())
  {
    return -fb.s1()*s1()-fb.s2()*s2()
           +fb.s3()*s3()+fb.s4()*s4();
  }
 
  template<typename ValueB>
  auto generalScalar(const LorentzSpinorBar<ValueB>& fb,
                                 Complex left, Complex right) const 
  -> decltype(left*fb.s1()*this->s1())
  {
    return  left*(fb.s1()*s1()+fb.s2()*s2())
         + right*(fb.s3()*s3()+fb.s4()*s4());
  }
 
  template<typename ValueB>
  auto sigma(const LorentzSpinorBar<ValueB>& fb) const 
    -> LorentzTensor<decltype(ValueB()*Value())> {
    typedef decltype(ValueB()*Value()) ResultT;
    LorentzTensor<ResultT> output;
    complex<ResultT> s11(fb.s1()*s1()),s22(fb.s2()*s2()),
      s33(fb.s3()*s3()),s44(fb.s4()*s4()),
      s12(fb.s1()*s2()),s21(fb.s2()*s1()),
      s34(fb.s3()*s4()),s43(fb.s4()*s3());
    Complex ii(0.,1.);
    ResultT zero;
    zero = ZERO;
    output.setTT(         zero         );
    output.setTX(-ii*( s12+s21-s34-s43));
    output.setTY(     -s12+s21+s34-s43 );
    output.setTZ(-ii*( s11-s22-s33+s44));
    output.setXT(      -output.tx()    );
    output.setXX(         zero         );
    output.setXY(      s11-s22+s33-s44 );
    output.setXZ(-ii*(-s12+s21-s34+s43));
    output.setYT(     -output.ty()     );
    output.setYX(     -output.xy()     );
    output.setYY(         zero         );
    output.setYZ(      s12+s21+s34+s43 );
    output.setZT(     -output.tz()     );
    output.setZX(     -output.xz()     );
    output.setZY(     -output.yz()     );
    output.setZZ(         zero         );
    return output;
  }
 
private:
  SpinorType _type;
 
  std::array<complex<Value>,4> _spin;
};
 
 
template <typename Value>
inline LorentzSpinor<double>
operator/(const LorentzSpinor<Value> & v, Value a) {
  return LorentzSpinor<double>(v.s1()/a, v.s2()/a, v.s3()/a, v.s4()/a,v.Type());
}
 
inline LorentzSpinor<double>
operator/(const LorentzSpinor<double> & v, Complex a) {
  return LorentzSpinor<double>(v.s1()/a, v.s2()/a, v.s3()/a, v.s4()/a,v.Type());
}
 
template <typename Value>
inline LorentzSpinor<Value> operator-(const LorentzSpinor<Value> & v) {
  return LorentzSpinor<Value>(-v.s1(),-v.s2(),-v.s3(),-v.s4(),v.Type());
}
 
template <typename ValueA, typename ValueB>
inline LorentzSpinor<ValueA>
operator+(LorentzSpinor<ValueA> a, const LorentzSpinor<ValueB> & b) {
  return a += b;
}
 
template <typename ValueA, typename ValueB>
inline LorentzSpinor<ValueA>
operator-(LorentzSpinor<ValueA> a, const LorentzSpinor<ValueB> & b) {
  return a -= b;
}
 
template <typename Value>
inline LorentzSpinor<Value>
operator*(const LorentzSpinor<Value> & a, double b) {
  return LorentzSpinor<Value>(a.s1()*b, a.s2()*b, a.s3()*b, a.s4()*b,a.Type());
}
 
template <typename Value>
inline LorentzSpinor<Value>
operator*(double b, LorentzSpinor<Value> a) {
  return a *= b;
}
  
template <typename Value>
inline LorentzSpinor<Value>
operator*(const LorentzSpinor<Value> & a, Complex b) {
  return LorentzSpinor<Value>(a.s1()*b, a.s2()*b, a.s3()*b, a.s4()*b,a.Type());
}
 
template <typename ValueA, typename ValueB>
inline auto operator*(complex<ValueB> a, const LorentzSpinor<ValueA> & v) 
  -> LorentzSpinor<decltype(a.real()*v.s1().real())>
{
  return {a*v.s1(), a*v.s2(), a*v.s3(), a*v.s4(),v.Type()};
}
 
template <typename ValueA, typename ValueB>
inline auto operator*(const LorentzSpinor<ValueA> & v, complex<ValueB> b) 
  -> LorentzSpinor<decltype(b.real()*v.s1().real())>
{
  return b*v;
}
 
template <typename ValueA, typename ValueB>
inline auto operator/(const LorentzSpinor<ValueA> & v, complex<ValueB> b) 
  -> LorentzSpinor<decltype(v.s1().real()/b.real())>
{
  return {v.s1()/b, v.s2()/b, v.s3()/b, v.s4()/b,v.Type()};
}
  
template <typename ValueA, typename ValueB>
inline auto operator*(ValueB a, const LorentzSpinor<ValueA> & v) 
  -> LorentzSpinor<decltype(a*v.s1().real())>
{
  return {a*v.s1(), a*v.s2(), a*v.s3(), a*v.s4(),v.Type()};
}
 
template <typename ValueA, typename ValueB>
inline auto operator*(const LorentzSpinor<ValueA> & v, ValueB b) 
  -> LorentzSpinor<decltype(b*v.s1().real())>
{
  return b*v;
}
 
template <typename ValueA, typename ValueB>
inline auto operator/(const LorentzSpinor<ValueA> & v, ValueB b) 
  -> LorentzSpinor<decltype(v.s1().real()/b)>
{
  return {v.s1()/b, v.s2()/b, v.s3()/b, v.s4()/b,v.Type()};
}
 
}
}
 
#ifndef ThePEG_TEMPLATES_IN_CC_FILE
#include "LorentzSpinor.tcc"
#endif 
 
#endif