diff --git a/Helicity/LorentzRank3Tensor.h b/Helicity/LorentzRank3Tensor.h
--- a/Helicity/LorentzRank3Tensor.h
+++ b/Helicity/LorentzRank3Tensor.h
@@ -1,278 +1,278 @@
 // -*- C++ -*-
 //
 // LorentzRank3Tensor.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_LorentzRank3Tensor_H
 #define ThePEG_LorentzRank3Tensor_H
 // This is the declaration of the LorentzRank3Tensor class.
 
 #include "ThePEG/Config/PhysicalQtyComplex.h"
 #include "ThePEG/Config/ThePEG.h"
 #include "LorentzTensor.h"
 
 namespace ThePEG {
 namespace Helicity {
 
 // compiler magic needs these pre-declarations to make friend templates work
 template<typename Value> class LorentzRank3Tensor;
 
 /**
  *  The LorentzRank3Tensor class is designed to implement the storage of a
  *  complex tensor to be used to representation the wavefunction of a
  *  spin-2 particle.
  *
  *  At the moment it only implements the storage of the tensor
  *  components but it is envisaged that it will be extended to include
  *  boost methods etc.
  *
  * @author Peter Richardson
  *
  */
 
 template<typename Value> 
 class LorentzRank3Tensor {
 
 public:
 
   /**
    * Default zero constructor.
    */
   LorentzRank3Tensor() = default;
 
   /**
    * Get components by indices.
    */
   complex<Value> operator () (int i, int j, int k) const {
     assert( i>=0 && i<=3 && j>=0 && j<=3 && k>=0 && k<=3);
     return _tensor[i][j][k];
   }
 
   /**
    * Set components by indices.
    */
   complex<Value> & operator () (int i, int j, int k) {
     assert( i>=0 && i<=3 && j>=0 && j<=3 && k>=0 && k<=3);
     return _tensor[i][j][k];
   }
   //@}
 
   /** @name Transformations. */
   //@{
   /**
    * Standard Lorentz boost specifying the components of the beta vector.
    */
   LorentzRank3Tensor & boost(double,double,double);
 
   /**
    * Standard Lorentz boost specifying the beta vector.
    */
   LorentzRank3Tensor<Value> & boost(const Boost & b) {
     return boost(b.x(), b.y(), b.z());
   }
 
   /**
    * General Lorentz transformation
    */
   LorentzRank3Tensor & transform(const SpinOneLorentzRotation & r){
     unsigned int ix,iy,iz,ixa,iya,iza;
     LorentzRank3Tensor<Value> output;
     complex<Value> temp;
     for(ix=0;ix<4;++ix) {
       for(iy=0;iy<4;++iy) {
 	for(iz=0;iz<4;++iz) {
 	  output(ix,iy,iz) = complex<Value>();
 	  for(ixa=0;ixa<4;++ixa) {
 	    for(iya=0;iya<4;++iya) {
 	      for(iza=0;iza<4;++iza)
 		 output(ix,iy,iz) += r(ix,ixa)*r(iy,iya)*r(iz,iza)*(*this)(ixa,iya,iza);
 	    }
 	  }
 	}
       }
     }
     *this=output;
     return *this;
   }
   
   /**
    * Return the complex conjugate.
    */
   LorentzRank3Tensor<Value> conjugate() {
     LorentzRank3Tensor<Value> output;
     for(unsigned int ix=0;ix<4;++ix) {
       for(unsigned int iy=0;iy<4;++iy) {
 	for(unsigned int iz=0;iz<4;++iz) {
 	  output(ix,iy,iz) = conj(output(ix,iy,iz));
 	}
       }
     }
     return output;
   }
 
   //@}
 
   /** @name Arithmetic operators. */
   //@{
   /**
    * Scaling with a complex number
    */
   LorentzRank3Tensor<Value> operator*=(Complex a) {
     for(int ix=0;ix<4;++ix)
       for(int iy=0;iy<4;++iy)
 	for(int iz=0;iz<4;++iz) _tensor[ix][iy][iz]*=a;
     return *this;
   }
 
   /**
    * Scalar product with other tensor
    */
   template <typename T, typename U>
   friend auto
   operator*(const LorentzRank3Tensor<T> & t, const LorentzRank3Tensor<U> & u) -> decltype(t.xx()*u.xx());
     
   /**
    * Addition.
    */
   LorentzRank3Tensor<Value> operator+(const LorentzRank3Tensor<Value> & in) const {
     LorentzRank3Tensor<Value> output;
     for(int ix=0;ix<4;++ix)
       for(int iy=0;iy<4;++iy)
 	for(int iz=0;iz<4;++iz) output(ix,iy,iz) = _tensor[ix][iy][iz] + in(ix,iy,iz);
   }
   
   /**
    * Subtraction.
    */
   LorentzRank3Tensor<Value> operator-(const LorentzRank3Tensor<Value> & in) const {
     LorentzRank3Tensor<Value> output;
     for(int ix=0;ix<4;++ix)
       for(int iy=0;iy<4;++iy)
 	for(int iz=0;iz<4;++iz) output(ix,iy,iz) = _tensor[ix][iy][iz] - in(ix,iy,iz);
   }
   
   /**
    * Dot product with the ith index
    */
   template<typename ValueB>
   auto dot(const LorentzVector<complex<ValueB> > & vec, unsigned int iloc) const 
-    -> LorentzTensor<decltype(vec.x()*Value())> {
+    -> LorentzTensor<decltype(ValueB()*Value())> {
     LorentzTensor<decltype(ValueB()*Value())> output;
     if(iloc==0) {
       for(unsigned int iy=0;iy<4;++iy) {
 	for(unsigned int iz=0;iz<4;++iz) {
 	  output(iy,iz) =
 	    vec.t()*_tensor[3][iy][iz] - vec.x()*_tensor[0][iy][iz] -
 	    vec.y()*_tensor[1][iy][iz] - vec.z()*_tensor[2][iy][iz];
 	}
       }
     }
     else if(iloc==1) {
       for(unsigned int iy=0;iy<4;++iy) {
 	for(unsigned int iz=0;iz<4;++iz) {
 	  output(iy,iz) =
 	    vec.t()*_tensor[iy][3][iz] - vec.x()*_tensor[iy][0][iz] -
 	    vec.y()*_tensor[iy][1][iz] - vec.z()*_tensor[iy][2][iz];
 	}
       }
     }
     else if(iloc==2) {
       for(unsigned int iy=0;iy<4;++iy) {
 	for(unsigned int iz=0;iz<4;++iz) {
 	  output(iy,iz) =
 	    vec.t()*_tensor[iy][iz][3] - vec.x()*_tensor[iy][iz][0] -
 	    vec.y()*_tensor[iy][iz][1] - vec.z()*_tensor[iy][iz][2];
 	}
       }
     }
     else
       assert(false);
     return output;
   }
   
   /**
    *  dot product with momentum
    */ 
   auto dot (const Lorentz5Momentum & vec,unsigned int iloc) const 
   -> LorentzTensor<decltype(vec.x()*Value())>
   {
     LorentzTensor<decltype(vec.x()*Value())> output;
     if(iloc==0) {
       for(unsigned int iy=0;iy<4;++iy) {
 	for(unsigned int iz=0;iz<4;++iz) {
 	  output(iy,iz) =
 	    vec.t()*_tensor[3][iy][iz] - vec.x()*_tensor[0][iy][iz] -
 	    vec.y()*_tensor[1][iy][iz] - vec.z()*_tensor[2][iy][iz];
 	}
       }
     }
     else if(iloc==1) {
       for(unsigned int iy=0;iy<4;++iy) {
 	for(unsigned int iz=0;iz<4;++iz) {
 	  output(iy,iz) =
 	    vec.t()*_tensor[iy][3][iz] - vec.x()*_tensor[iy][0][iz] -
 	    vec.y()*_tensor[iy][1][iz] - vec.z()*_tensor[iy][2][iz];
 	}
       }
     }
     else if(iloc==2) {
       for(unsigned int iy=0;iy<4;++iy) {
 	for(unsigned int iz=0;iz<4;++iz) {
 	  output(iy,iz) =
 	    vec.t()*_tensor[iy][iz][3] - vec.x()*_tensor[iy][iz][0] -
 	    vec.y()*_tensor[iy][iz][1] - vec.z()*_tensor[iy][iz][2];
 	}
       }
     }
     else
       assert(false);
     return output;
   } 
   //@}
 
 private:
 
   /**
    * The components.
    */
   std::array<std::array<std::array<complex<Value>,4>,4>,4> _tensor;
 
 };
 
 /**
  * Multiplication by a complex number.
  */
 template<typename T, typename U> 
 inline auto
 operator*(complex<U> a, const LorentzRank3Tensor<T> & t) -> LorentzRank3Tensor<decltype(a.real()*t.xx().real())> {
   LorentzRank3Tensor<decltype(a.real()*t.xx().real())> output;
   for(int ix=0;ix<4;++ix)
     for(int iy=0;iy<4;++iy)
       for(int iz=0;iz<4;++iz) output(ix,iy,iz) = a*t(ix,iy,iz);
   return output;
 }
 
 /**
  * Multiplication by a complex number.
  */
 template<typename T, typename U> 
 inline auto
 operator*(const LorentzRank3Tensor<T> & t,complex<U> a) -> LorentzRank3Tensor<decltype(a.real()*t.xx().real())> {
   LorentzRank3Tensor<decltype(a.real()*t.xx().real())> output;
   for(int ix=0;ix<4;++ix)
     for(int iy=0;iy<4;++iy)
       for(int iz=0;iz<4;++iz) output(ix,iy,iz) = a*t(ix,iy,iz);
   return output;
 }
 
 }
 }
 
 #ifndef ThePEG_TEMPLATES_IN_CC_FILE
 #include "LorentzRank3Tensor.tcc"
 #endif 
 
 #endif
diff --git a/Helicity/LorentzTensor.h b/Helicity/LorentzTensor.h
--- a/Helicity/LorentzTensor.h
+++ b/Helicity/LorentzTensor.h
@@ -1,572 +1,572 @@
 // -*- 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;
 
 /**
  *  The LorentzTensor class is designed to implement the storage of a
  *  complex tensor to be used to representation the wavefunction of a
  *  spin-2 particle.
  *
  *  At the moment it only implements the storage of the tensor
  *  components but it is envisaged that it will be extended to include
  *  boost methods etc.
  *
  * @author Peter Richardson
  *
  */
 
 template<typename Value> 
 class LorentzTensor {
 
 public:
 
   /** @name Standard constructors and destructors. */
   //@{
   /**
    * Default zero constructor.
    */
   LorentzTensor() = default;
 
   /**
    * Constructor specifyign all components.
    */
   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}} }} {}
   /**
    * Constructor in terms of two polarization vectors.
    */
   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());
   }
   //@}
 
   /** @name Access individual components. */
   //@{
   /**
    * Get x,x component.
    */
   complex<Value> xx() const {return _tensor[0][0];}
 
   /**
    * Get y,x component.
    */
   complex<Value> yx() const {return _tensor[1][0];}
   /**
    * Get z,x component.
    */
   complex<Value> zx() const {return _tensor[2][0];}
 
   /**
    * Get t,x component.
    */
   complex<Value> tx() const {return _tensor[3][0];}
 
   /**
    * Get x,y component.
    */
   complex<Value> xy() const {return _tensor[0][1];}
 
   /**
    * Get y,y component.
    */
   complex<Value> yy() const {return _tensor[1][1];}
 
   /**
    * Get z,y component.
    */
   complex<Value> zy() const {return _tensor[2][1];}
 
   /**
    * Get t,y component.
    */
   complex<Value> ty() const {return _tensor[3][1];}
 
   /**
    * Get x,z component.
    */
   complex<Value> xz() const {return _tensor[0][2];}
 
   /**
    * Get y,z component.
    */
   complex<Value> yz() const {return _tensor[1][2];}
 
   /**
    * Get z,z component.
    */
   complex<Value> zz() const {return _tensor[2][2];}
 
   /**
    * Get t,z component.
    */
   complex<Value> tz() const {return _tensor[3][2];}
 
   /**
    * Get x,t component.
    */
   complex<Value> xt() const {return _tensor[0][3];}
 
   /**
    * Get y,t component.
    */
   complex<Value> yt() const {return _tensor[1][3];}
 
   /**
    * Get z,t component.
    */
   complex<Value> zt() const {return _tensor[2][3];}
 
   /**
    * Get t,t component.
    */
   complex<Value> tt() const {return _tensor[3][3];}
 
   /**
    * Set x,x component.
    */
   void setXX(complex<Value> a) {_tensor[0][0]=a;}
 
   /**
    * Set y,x component.
    */
   void setYX(complex<Value> a) {_tensor[1][0]=a;}
 
   /**
    * Set z,x component.
    */
   void setZX(complex<Value> a) {_tensor[2][0]=a;}
 
   /**
    * Set t,x component.
    */
   void setTX(complex<Value> a) {_tensor[3][0]=a;}
 
   /**
    * Set x,y component.
    */
   void setXY(complex<Value> a) {_tensor[0][1]=a;}
 
   /**
    * Set y,y component.
    */
   void setYY(complex<Value> a) {_tensor[1][1]=a;}
 
   /**
    * Set z,y component.
    */
   void setZY(complex<Value> a) {_tensor[2][1]=a;}
 
   /**
    * Set t,y component.
    */
   void setTY(complex<Value> a) {_tensor[3][1]=a;}
 
   /**
    * Set x,z component.
    */
   void setXZ(complex<Value> a) {_tensor[0][2]=a;}
 
   /**
    * Set y,z component.
    */
   void setYZ(complex<Value> a) {_tensor[1][2]=a;}
 
   /**
    * Set z,z component.
    */
   void setZZ(complex<Value> a) {_tensor[2][2]=a;}
 
   /**
    * Set t,z component.
    */
   void setTZ(complex<Value> a) {_tensor[3][2]=a;}
 
   /**
    * Set x,t component.
    */
   void setXT(complex<Value> a) {_tensor[0][3]=a;}
 
   /**
    * Set y,t component.
    */
   void setYT(complex<Value> a) {_tensor[1][3]=a;}
 
   /**
    * Set z,t component.
    */
   void setZT(complex<Value> a) {_tensor[2][3]=a;}
 
   /**
    * Set t,t component.
    */
   void setTT(complex<Value> a) {_tensor[3][3]=a;}
 
   /**
    * Get components by indices.
    */
   complex<Value> operator () (int i, int j) const {
     assert( i>=0 && i<=3 && j>=0 && j<=3);
     return _tensor[i][j];
   }
 
   /**
    * Set components by indices.
    */
   complex<Value> & operator () (int i, int j) {
     assert( i>=0 && i<=3 && j>=0 && j<=3);
     return _tensor[i][j];
   }
   //@}
 
   /** @name Transformations. */
   //@{
   /**
    * Standard Lorentz boost specifying the components of the beta vector.
    */
   LorentzTensor & boost(double,double,double);
 
   /**
    * Standard Lorentz boost specifying the beta vector.
    */
   LorentzTensor<Value> & boost(const Boost & b) {
     return boost(b.x(), b.y(), b.z());
   }
 
   /**
    * General Lorentz transformation
    */
   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;
   }
   
   /**
    * Return the complex conjugate.
    */
   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()));
   }
 
   //@}
 
   /** @name Arithmetic operators. */
   //@{
   /**
    * Scaling with a complex number
    */
   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;
   }
 
   /**
    * Scalar product with other tensor
    */
   template <typename T, typename U>
   friend auto
   operator*(const LorentzTensor<T> & t, const LorentzTensor<U> & u)
   -> decltype(t.xx()*u.xx())
   ;
     
   /**
    * Addition.
    */
   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());
   }
   
   /**
    * Subtraction.
    */
   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());
   }
 
   /**
    * Trace
    */
   complex<Value> trace() const {
     return _tensor[3][3]-_tensor[0][0]-_tensor[1][1]-_tensor[2][2];
   }
 
   /**
    *  Inner product with another tensor
    */
   template<typename ValueB>
   auto innerProduct(const LorentzTensor<ValueB> & ten) const {
     LorentzTensor<decltype(ten.xx().real()*this->xx().real())> output;
-    for(unsigned int ix=0;ix<3;++ix) {
-      for(unsigned int iy=0;iy<3;++iy) {
+    for(unsigned int ix=0;ix<4;++ix) {
+      for(unsigned int iy=0;iy<4;++iy) {
 	output(ix,iy) = _tensor[ix][3]*ten(3,iy);
 	for(unsigned int iz=0;iz<3;++iz) {
 	  output(ix,iy) -= _tensor[ix][iz]*ten(iz,iy);
 	}
       }
     }
     return output;
   }
 
   /**
    *  Outer product with another tensor
    */
   template<typename ValueB>
   auto outerProduct(const LorentzTensor<ValueB> & ten) const {
     LorentzTensor<decltype(ten.xx().real()*this->xx().real())> output;
-    for(unsigned int ix=0;ix<3;++ix) {
-      for(unsigned int iy=0;iy<3;++iy) {
+    for(unsigned int ix=0;ix<4;++ix) {
+      for(unsigned int iy=0;iy<4;++iy) {
 	output(ix,iy) = _tensor[3][ix]*ten(iy,3);
 	for(unsigned int iz=0;iz<3;++iz) {
 	  output(ix,iy) -= _tensor[iz][ix]*ten(iy,iz);
 	}
       }
     }
     return output;
   }
 
   //@}
 
   /**
    *  Various dot products
    */
   //@{
   /**
    *  First index dot product with polarization vector
    */
   template<typename ValueB>
   auto preDot (const LorentzVector<complex<ValueB> > & 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;
   }
 
   /**
    *  Second index dot product with polarization vector
    */
   template<typename ValueB>
   auto postDot(const LorentzVector<complex<ValueB> > & 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;
   }
 
   /**
    *  First index dot product with momentum
    */ 
   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;
   }
 
   /**
    *  Second index dot product with momentum
    */ 
   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:
 
   /**
    * The components.
    */
   std::array<std::array<complex<Value>,4>,4> _tensor;
 
 };
 
 /**
  * Multiplication by a complex number.
  */
 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()};
 }
 
 /**
  * Multiplication by a complex number.
  */
 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()};
 }
 
 /**
  * Multiply a LorentzVector by a LorentzTensor.
  */
 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;
 }
 
 /**
  * Multiply a LorentzTensor by a LorentzVector.
  */
 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;
 }
 
 /**
  * Multiply a LorentzTensor by a LorentzTensor
  */
 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
diff --git a/Helicity/epsilon.h b/Helicity/epsilon.h
--- a/Helicity/epsilon.h
+++ b/Helicity/epsilon.h
@@ -1,165 +1,227 @@
 // -*- C++ -*-
 //
 // epsilon.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_epsilon_H
 #define ThePEG_epsilon_H
 //
 // This is the declaration of the epsilon class.
 
 #include "ThePEG/Vectors/LorentzVector.h"
 #include "LorentzTensor.h"
 
 namespace ThePEG {
 namespace Helicity {
 
 /** \ingroup Helicity
  *  \author Peter Richardson
  *
  *  This class is designed to combine 5-momenta and polarization 
  *  vectors together with the result being the product with the 
  *  eps function. The class is purely static and contains no data.
  *
  *  @see LorentzPolarizationVector
  *  @see Lorentz5Vector
  */
 
   /**
    *  Return the product 
    *  \f$\epsilon^{\alpha\beta\gamma\delta}v_{1\alpha}v_{2\beta}v_{3\gamma}v_{4\delta}\f$.
    * @param a The first  vector \f$v_{1\alpha}\f$.
    * @param b The second vector \f$v_{2\beta}\f$.
    * @param c The third  vector \f$v_{3\gamma}\f$.
    * @param d The fourth vector \f$v_{4\delta}\f$.
    * @return The product 
    * \f$\epsilon^{\alpha\beta\gamma\delta}v_{1\alpha}v_{2\beta}v_{3\gamma}v_{4\delta}\f$
    */
   template <typename A, typename B, typename C, typename D>
   auto epsilon(const LorentzVector<A> & a,
                const LorentzVector<B> & b,
                const LorentzVector<C> & c,
                const LorentzVector<D> & d) 
     -> decltype(a.x()*b.y()*c.z()*d.t())
   {
     auto diffxy = a.x() * b.y()  -  a.y() * b.x();
     auto diffxz = a.x() * b.z()  -  a.z() * b.x();
     auto diffxt = a.x() * b.t()  -  a.t() * b.x();
     auto diffyz = a.y() * b.z()  -  a.z() * b.y();
     auto diffyt = a.y() * b.t()  -  a.t() * b.y();
     auto diffzt = a.z() * b.t()  -  a.t() * b.z();
     
     auto diff2xy = c.x() * d.y()  -  c.y() * d.x();
     auto diff2xz = c.x() * d.z()  -  c.z() * d.x();
     auto diff2xt = c.x() * d.t()  -  c.t() * d.x();
     auto diff2yz = c.y() * d.z()  -  c.z() * d.y();
     auto diff2yt = c.y() * d.t()  -  c.t() * d.y();
     auto diff2zt = c.z() * d.t()  -  c.t() * d.z();
  
     return
       diff2yz*diffxt + diff2zt*diffxy - diff2yt*diffxz - 
       diff2xz*diffyt + diff2xt*diffyz + diff2xy*diffzt;
   }
   
   /**
    *  Return the product 
    *  \f$\epsilon^{\mu\alpha\beta\gamma}v_{1\alpha}v_{2\beta}v_{3\gamma}\f$.
    * @param a The first  vector \f$v_{1\alpha}\f$.
    * @param b The second vector \f$v_{2\beta}\f$.
    * @param c The third  vector \f$v_{3\gamma}\f$.
    * @return The product 
    * \f$\epsilon^{\mu\alpha\beta\gamma}v_{1\alpha}v_{2\beta}v_{3\gamma}\f$.
    */
   template <typename A, typename B, typename C>
   auto epsilon(const LorentzVector<A> & a,
                const LorentzVector<B> & b,
                const LorentzVector<C> & c) 
   -> LorentzVector<decltype(a.x()*b.y()*c.z())>
   {
     auto diffxy = a.x() * b.y()  -  a.y() * b.x();
     auto diffxz = a.x() * b.z()  -  a.z() * b.x();
     auto diffxt = a.x() * b.t()  -  a.t() * b.x();
     auto diffyz = a.y() * b.z()  -  a.z() * b.y();
     auto diffyt = a.y() * b.t()  -  a.t() * b.y();
     auto diffzt = a.z() * b.t()  -  a.t() * b.z();
 
     using ResultType = LorentzVector<decltype(a.x()*b.x()*c.x())>;    
     ResultType result;
     result.setX( c.z() * diffyt  - c.t() * diffyz  - c.y() * diffzt); 
     result.setY( c.t() * diffxz  - c.z() * diffxt  + c.x() * diffzt);
     result.setZ(-c.t() * diffxy  + c.y() * diffxt  - c.x() * diffyt);
     result.setT(-c.z() * diffxy  + c.y() * diffxz  - c.x() * diffyz);
     
     return result;
   }
   
   /**
    *  Return the product 
    *  \f$\epsilon^{\mu\nu\alpha\beta}v_{1\alpha}v_{2\beta}\f$.
    * @param a The first  vector \f$v_{1\alpha}\f$.
    * @param b The second vector \f$v_{2\beta}\f$.
    * @return The product
    * \f$\epsilon^{\mu\nu\alpha\beta\gamma}v_{1\alpha}v_{2\beta}\f$.
    */
   template <typename A, typename B>
   auto epsilon(const LorentzVector<complex<A> > & a,
                const LorentzVector<complex<B> > & b)
     -> LorentzTensor<decltype(a.x().real()*b.y().real())>
   {
     auto diffxy = a.x() * b.y()  -  a.y() * b.x();
     auto diffxz = a.x() * b.z()  -  a.z() * b.x();
     auto diffxt = a.x() * b.t()  -  a.t() * b.x();
     auto diffyz = a.y() * b.z()  -  a.z() * b.y();
     auto diffyt = a.y() * b.t()  -  a.t() * b.y();
     auto diffzt = a.z() * b.t()  -  a.t() * b.z();
     complex<decltype(a.x().real()*b.x().real())> zero(ZERO);
 
     using ResultType = LorentzTensor<decltype(a.x().real()*b.x().real())>;
     ResultType result;
     result.setTT(  zero ); result.setTX( diffyz); result.setTY(-diffxz); result.setTZ( diffxy);
     result.setXT(-diffyz); result.setXX(  zero ); result.setXY( diffzt); result.setXZ(-diffyt);
     result.setYT( diffxz); result.setYX(-diffzt); result.setYY(  zero ); result.setYZ( diffxt);
     result.setZT(-diffxy); result.setZX( diffyt); result.setZY(-diffxt); result.setZZ(  zero );
 
     return result;
   }
+  
+  /**
+   *  Return the product 
+   *  \f$\epsilon^{\mu\nu\alpha\beta}v_{1\alpha}v_{2\beta}\f$.
+   * @param a The first  vector \f$v_{1\alpha}\f$.
+   * @param b The second vector \f$v_{2\beta}\f$.
+   * @return The product
+   * \f$\epsilon^{\mu\nu\alpha\beta\gamma}v_{1\alpha}v_{2\beta}\f$.
+   */
+  template <typename A, typename B>
+  auto epsilon(const LorentzVector<A>           & a,
+               const LorentzVector<complex<B> > & b)
+    -> LorentzTensor<decltype(a.x()*b.y().real())>
+  {
+    auto diffxy = a.x() * b.y()  -  a.y() * b.x();
+    auto diffxz = a.x() * b.z()  -  a.z() * b.x();
+    auto diffxt = a.x() * b.t()  -  a.t() * b.x();
+    auto diffyz = a.y() * b.z()  -  a.z() * b.y();
+    auto diffyt = a.y() * b.t()  -  a.t() * b.y();
+    auto diffzt = a.z() * b.t()  -  a.t() * b.z();
+    complex<decltype(a.x()*b.x().real())> zero(ZERO);
+
+    using ResultType = LorentzTensor<decltype(a.x()*b.x().real())>;
+    ResultType result;
+    result.setTT(  zero ); result.setTX( diffyz); result.setTY(-diffxz); result.setTZ( diffxy);
+    result.setXT(-diffyz); result.setXX(  zero ); result.setXY( diffzt); result.setXZ(-diffyt);
+    result.setYT( diffxz); result.setYX(-diffzt); result.setYY(  zero ); result.setYZ( diffxt);
+    result.setZT(-diffxy); result.setZX( diffyt); result.setZY(-diffxt); result.setZZ(  zero );
+
+    return result;
+  }
+
+  /**
+   *  Return the product 
+   *  \f$\epsilon^{\mu\nu\alpha\beta}v_{1\alpha}v_{2\beta}\f$.
+   * @param a The first  vector \f$v_{1\alpha}\f$.
+   * @param b The second vector \f$v_{2\beta}\f$.
+   * @return The product
+   * \f$\epsilon^{\mu\nu\alpha\beta\gamma}v_{1\alpha}v_{2\beta}\f$.
+   */
+  template <typename A, typename B>
+  auto epsilon(const LorentzVector<complex<A> > & a,
+               const LorentzVector<B>           & b)
+    -> LorentzTensor<decltype(a.x().real()*b.y())>
+  {
+    auto diffxy = a.x() * b.y()  -  a.y() * b.x();
+    auto diffxz = a.x() * b.z()  -  a.z() * b.x();
+    auto diffxt = a.x() * b.t()  -  a.t() * b.x();
+    auto diffyz = a.y() * b.z()  -  a.z() * b.y();
+    auto diffyt = a.y() * b.t()  -  a.t() * b.y();
+    auto diffzt = a.z() * b.t()  -  a.t() * b.z();
+    complex<decltype(a.x().real()*b.x())> zero(ZERO);
+
+    using ResultType = LorentzTensor<decltype(a.x().real()*b.x())>;
+    ResultType result;
+    result.setTT(  zero ); result.setTX( diffyz); result.setTY(-diffxz); result.setTZ( diffxy);
+    result.setXT(-diffyz); result.setXX(  zero ); result.setXY( diffzt); result.setXZ(-diffyt);
+    result.setYT( diffxz); result.setYX(-diffzt); result.setYY(  zero ); result.setYZ( diffxt);
+    result.setZT(-diffxy); result.setZX( diffyt); result.setZY(-diffxt); result.setZZ(  zero );
+
+    return result;
+  }
 
   /**
    *  Return the product
    *  \f$\epsilon^{\mu\nu\alpha\beta}v_{1\alpha}v_{2\beta}\f$.
    * @param a The first  vector \f$v_{1\alpha}\f$.
    * @param b The second vector \f$v_{2\beta}\f$.
    * @return The product
    * \f$\epsilon^{\mu\nu\alpha\beta\gamma}v_{1\alpha}v_{2\beta}\f$.
    */
   template <typename A, typename B>
   auto epsilon(const LorentzVector<A> & a,
                const LorentzVector<B> & b)
     -> LorentzTensor<decltype(a.x()*b.y())>
   {
     auto diffxy = a.x() * b.y()  -  a.y() * b.x();
     auto diffxz = a.x() * b.z()  -  a.z() * b.x();
     auto diffxt = a.x() * b.t()  -  a.t() * b.x();
     auto diffyz = a.y() * b.z()  -  a.z() * b.y();
     auto diffyt = a.y() * b.t()  -  a.t() * b.y();
     auto diffzt = a.z() * b.t()  -  a.t() * b.z();
     complex<decltype(a.x()*b.x())> zero(ZERO);
 
     using ResultType = LorentzTensor<decltype(a.x()*b.x())>;
     ResultType result;
     result.setTT(  zero ); result.setTX( diffyz); result.setTY(-diffxz); result.setTZ( diffxy);
     result.setXT(-diffyz); result.setXX(  zero ); result.setXY( diffzt); result.setXZ(-diffyt);
     result.setYT( diffxz); result.setYX(-diffzt); result.setYY(  zero ); result.setYZ( diffxt);
     result.setZT(-diffxy); result.setZX( diffyt); result.setZY(-diffxt); result.setZZ(  zero );
 
     return result;
   }
 
 
 }
 }
 
 #endif /* ThePEG_epsilon_H */