ThePEG/Helicity/LorentzTensor.h

0001 // -*- C++ -*-
0002 //
0003 // LorentzTensor.h is a part of ThePEG - Toolkit for HEP Event Generation
0004 // Copyright (C) 2003-2019 Peter Richardson, Leif Lonnblad
0005 //
0006 // ThePEG is licenced under version 3 of the GPL, see COPYING for details.
0007 // Please respect the MCnet academic guidelines, see GUIDELINES for details.
0008 //
0009 #ifndef ThePEG_LorentzTensor_H
0010 #define ThePEG_LorentzTensor_H
0011 // This is the declaration of the LorentzTensor class.
0012
0013 #include "ThePEG/Config/PhysicalQtyComplex.h"
0014 #include "ThePEG/Config/ThePEG.h"
0015 #include "LorentzPolarizationVector.h"
0016
0017 namespace ThePEG {
0018 namespace Helicity {
0019
0020 // compiler magic needs these pre-declarations to make friend templates work
0021 template<typename Value> class LorentzTensor;
0022
0023 /**
0024  *  The LorentzTensor class is designed to implement the storage of a
0025  *  complex tensor to be used to representation the wavefunction of a
0026  *  spin-2 particle.
0027  *
0028  *  At the moment it only implements the storage of the tensor
0029  *  components but it is envisaged that it will be extended to include
0030  *  boost methods etc.
0031  *
0032  * @author Peter Richardson
0033  *
0034  */
0035
0036 template<typename Value>
0037 class LorentzTensor {
0038
0039 public:
0040
0041   /** @name Standard constructors and destructors. */
0042   //@{
0043   /**
0044    * Default zero constructor.
0045    */
0046   LorentzTensor() = default;
0047
0048   /**
0049    * Constructor specifyign all components.
0050    */
0051   LorentzTensor(complex<Value> xx, complex<Value> xy,
0052                 complex<Value> xz, complex<Value> xt,
0053                 complex<Value> yx, complex<Value> yy,
0054                 complex<Value> yz, complex<Value> yt,
0055                 complex<Value> zx, complex<Value> zy,
0056                 complex<Value> zz, complex<Value> zt,
0057                 complex<Value> tx, complex<Value> ty,
0058                 complex<Value> tz, complex<Value> tt)
0059   : _tensor{{ {{xx,xy,xz,xt}},
0060               {{yx,yy,yz,yt}},
0061               {{zx,zy,zz,zt}},
0062               {{tx,ty,tz,tt}} }} {}
0063   /**
0064    * Constructor in terms of two polarization vectors.
0065    */
0066   LorentzTensor(const LorentzPolarizationVector & p,
0067                 const LorentzPolarizationVector & q) {
0068     setXX(p.x() * q.x()); setYX(p.y() * q.x());
0069     setZX(p.z() * q.x()); setTX(p.t() * q.x());
0070     setXY(p.x() * q.y()); setYY(p.y() * q.y());
0071     setZY(p.z() * q.y()); setTY(p.t() * q.y());
0072     setXZ(p.x() * q.z()); setYZ(p.y() * q.z());
0073     setZZ(p.z() * q.z()); setTZ(p.t() * q.z());
0074     setXT(p.x() * q.t()); setYT(p.y() * q.t());
0075     setZT(p.z() * q.t()); setTT(p.t() * q.t());
0076   }
0077   //@}
0078
0079   /** @name Access individual components. */
0080   //@{
0081   /**
0082    * Get x,x component.
0083    */
0084   complex<Value> xx() const {return _tensor[0][0];}
0085
0086   /**
0087    * Get y,x component.
0088    */
0089   complex<Value> yx() const {return _tensor[1][0];}
0090   /**
0091    * Get z,x component.
0092    */
0093   complex<Value> zx() const {return _tensor[2][0];}
0094
0095   /**
0096    * Get t,x component.
0097    */
0098   complex<Value> tx() const {return _tensor[3][0];}
0099
0100   /**
0101    * Get x,y component.
0102    */
0103   complex<Value> xy() const {return _tensor[0][1];}
0104
0105   /**
0106    * Get y,y component.
0107    */
0108   complex<Value> yy() const {return _tensor[1][1];}
0109
0110   /**
0111    * Get z,y component.
0112    */
0113   complex<Value> zy() const {return _tensor[2][1];}
0114
0115   /**
0116    * Get t,y component.
0117    */
0118   complex<Value> ty() const {return _tensor[3][1];}
0119
0120   /**
0121    * Get x,z component.
0122    */
0123   complex<Value> xz() const {return _tensor[0][2];}
0124
0125   /**
0126    * Get y,z component.
0127    */
0128   complex<Value> yz() const {return _tensor[1][2];}
0129
0130   /**
0131    * Get z,z component.
0132    */
0133   complex<Value> zz() const {return _tensor[2][2];}
0134
0135   /**
0136    * Get t,z component.
0137    */
0138   complex<Value> tz() const {return _tensor[3][2];}
0139
0140   /**
0141    * Get x,t component.
0142    */
0143   complex<Value> xt() const {return _tensor[0][3];}
0144
0145   /**
0146    * Get y,t component.
0147    */
0148   complex<Value> yt() const {return _tensor[1][3];}
0149
0150   /**
0151    * Get z,t component.
0152    */
0153   complex<Value> zt() const {return _tensor[2][3];}
0154
0155   /**
0156    * Get t,t component.
0157    */
0158   complex<Value> tt() const {return _tensor[3][3];}
0159
0160   /**
0161    * Set x,x component.
0162    */
0163   void setXX(complex<Value> a) {_tensor[0][0]=a;}
0164
0165   /**
0166    * Set y,x component.
0167    */
0168   void setYX(complex<Value> a) {_tensor[1][0]=a;}
0169
0170   /**
0171    * Set z,x component.
0172    */
0173   void setZX(complex<Value> a) {_tensor[2][0]=a;}
0174
0175   /**
0176    * Set t,x component.
0177    */
0178   void setTX(complex<Value> a) {_tensor[3][0]=a;}
0179
0180   /**
0181    * Set x,y component.
0182    */
0183   void setXY(complex<Value> a) {_tensor[0][1]=a;}
0184
0185   /**
0186    * Set y,y component.
0187    */
0188   void setYY(complex<Value> a) {_tensor[1][1]=a;}
0189
0190   /**
0191    * Set z,y component.
0192    */
0193   void setZY(complex<Value> a) {_tensor[2][1]=a;}
0194
0195   /**
0196    * Set t,y component.
0197    */
0198   void setTY(complex<Value> a) {_tensor[3][1]=a;}
0199
0200   /**
0201    * Set x,z component.
0202    */
0203   void setXZ(complex<Value> a) {_tensor[0][2]=a;}
0204
0205   /**
0206    * Set y,z component.
0207    */
0208   void setYZ(complex<Value> a) {_tensor[1][2]=a;}
0209
0210   /**
0211    * Set z,z component.
0212    */
0213   void setZZ(complex<Value> a) {_tensor[2][2]=a;}
0214
0215   /**
0216    * Set t,z component.
0217    */
0218   void setTZ(complex<Value> a) {_tensor[3][2]=a;}
0219
0220   /**
0221    * Set x,t component.
0222    */
0223   void setXT(complex<Value> a) {_tensor[0][3]=a;}
0224
0225   /**
0226    * Set y,t component.
0227    */
0228   void setYT(complex<Value> a) {_tensor[1][3]=a;}
0229
0230   /**
0231    * Set z,t component.
0232    */
0233   void setZT(complex<Value> a) {_tensor[2][3]=a;}
0234
0235   /**
0236    * Set t,t component.
0237    */
0238   void setTT(complex<Value> a) {_tensor[3][3]=a;}
0239
0240   /**
0241    * Get components by indices.
0242    */
0243   complex<Value> operator () (int i, int j) const {
0244     assert( i>=0 && i<=3 && j>=0 && j<=3);
0245     return _tensor[i][j];
0246   }
0247
0248   /**
0249    * Set components by indices.
0250    */
0251   complex<Value> & operator () (int i, int j) {
0252     assert( i>=0 && i<=3 && j>=0 && j<=3);
0253     return _tensor[i][j];
0254   }
0255   //@}
0256
0257   /** @name Transformations. */
0258   //@{
0259   /**
0260    * Standard Lorentz boost specifying the components of the beta vector.
0261    */
0262   LorentzTensor & boost(double,double,double);
0263
0264   /**
0265    * Standard Lorentz boost specifying the beta vector.
0266    */
0267   LorentzTensor<Value> & boost(const Boost & b) {
0268     return boost(b.x(), b.y(), b.z());
0269   }
0270
0271   /**
0272    * General Lorentz transformation
0273    */
0274   LorentzTensor & transform(const SpinOneLorentzRotation & r){
0275     unsigned int ix,iy,ixa,iya;
0276     LorentzTensor<Value> output;
0277     complex<Value> temp;
0278     for(ix=0;ix<4;++ix) {
0279       for(iy=0;iy<4;++iy) {
0280         temp=complex<Value>();
0281         for(ixa=0;ixa<4;++ixa) {
0282           for(iya=0;iya<4;++iya)
0283             temp+=r(ix,ixa)*r(iy,iya)*(*this)(ixa,iya);
0284         }
0285         output(ix,iy)=temp;
0286       }
0287     }
0288     *this=output;
0289     return *this;
0290   }
0291
0292   /**
0293    * Return the complex conjugate.
0294    */
0295   LorentzTensor<Value> conjugate() {
0296     return LorentzTensor<Value>(conj(xx()), conj(xy()), conj(xz()), conj(xt()),
0297                                 conj(yx()), conj(yy()), conj(yz()), conj(yt()),
0298                                 conj(zx()), conj(zy()), conj(zz()), conj(zt()),
0299                                 conj(tx()), conj(ty()), conj(tz()), conj(tt()));
0300   }
0301
0302   //@}
0303
0304   /** @name Arithmetic operators. */
0305   //@{
0306   /**
0307    * Scaling with a complex number
0308    */
0309   LorentzTensor<Value> operator*=(Complex a) {
0310     for(int ix=0;ix<4;++ix)
0311       for(int iy=0;iy<4;++iy) _tensor[ix][iy]*=a;
0312     return *this;
0313   }
0314
0315   /**
0316    * Scalar product with other tensor
0317    */
0318   template <typename T, typename U>
0319   friend auto
0320   operator*(const LorentzTensor<T> & t, const LorentzTensor<U> & u)
0321   -> decltype(t.xx()*u.xx())
0322   ;
0323
0324   /**
0325    * Addition.
0326    */
0327   LorentzTensor<Value> operator+(const LorentzTensor<Value> & in) const {
0328     return LorentzTensor<Value>(xx()+in.xx(),xy()+in.xy(),xz()+in.xz(),xt()+in.xt(),
0329                                 yx()+in.yx(),yy()+in.yy(),yz()+in.yz(),yt()+in.yt(),
0330                                 zx()+in.zx(),zy()+in.zy(),zz()+in.zz(),zt()+in.zt(),
0331                                 tx()+in.tx(),ty()+in.ty(),tz()+in.tz(),tt()+in.tt());
0332   }
0333
0334   /**
0335    * Subtraction.
0336    */
0337   LorentzTensor<Value> operator-(const LorentzTensor<Value> & in) const {
0338     return LorentzTensor<Value>(xx()-in.xx(),xy()-in.xy(),xz()-in.xz(),xt()-in.xt(),
0339                                 yx()-in.yx(),yy()-in.yy(),yz()-in.yz(),yt()-in.yt(),
0340                                 zx()-in.zx(),zy()-in.zy(),zz()-in.zz(),zt()-in.zt(),
0341                                 tx()-in.tx(),ty()-in.ty(),tz()-in.tz(),tt()-in.tt());
0342   }
0343
0344   /**
0345    * Trace
0346    */
0347   complex<Value> trace() const {
0348     return _tensor[3][3]-_tensor[0][0]-_tensor[1][1]-_tensor[2][2];
0349   }
0350
0351   /**
0352    *  Inner product with another tensor
0353    */
0354   template<typename ValueB>
0355   auto innerProduct(const LorentzTensor<ValueB> & ten) const
0356     -> LorentzTensor<decltype(ten.xx().real()*this->xx().real())> {
0357     LorentzTensor<decltype(ten.xx().real()*this->xx().real())> output;
0358     for(unsigned int ix=0;ix<4;++ix) {
0359       for(unsigned int iy=0;iy<4;++iy) {
0360     output(ix,iy) = _tensor[ix][3]*ten(3,iy);
0361     for(unsigned int iz=0;iz<3;++iz) {
0362       output(ix,iy) -= _tensor[ix][iz]*ten(iz,iy);
0363     }
0364       }
0365     }
0366     return output;
0367   }
0368
0369   /**
0370    *  Outer product with another tensor
0371    */
0372   template<typename ValueB>
0373   auto outerProduct(const LorentzTensor<ValueB> & ten) const
0374     -> LorentzTensor<decltype(ten.xx().real()*this->xx().real())> {
0375     LorentzTensor<decltype(ten.xx().real()*this->xx().real())> output;
0376     for(unsigned int ix=0;ix<4;++ix) {
0377       for(unsigned int iy=0;iy<4;++iy) {
0378     output(ix,iy) = _tensor[3][ix]*ten(iy,3);
0379     for(unsigned int iz=0;iz<3;++iz) {
0380       output(ix,iy) -= _tensor[iz][ix]*ten(iy,iz);
0381     }
0382       }
0383     }
0384     return output;
0385   }
0386
0387   //@}
0388
0389   /**
0390    *  Various dot products
0391    */
0392   //@{
0393   /**
0394    *  First index dot product with polarization vector
0395    */
0396   template<typename ValueB>
0397   auto preDot (const LorentzVector<complex<ValueB> > & vec) const
0398     -> LorentzVector<decltype(vec.x()*this->xx())> {
0399     LorentzVector<decltype(vec.x()*this->xx())> output;
0400     output.setX(vec.t()*_tensor[3][0]-vec.x()*_tensor[0][0]-
0401                 vec.y()*_tensor[1][0]-vec.z()*_tensor[2][0]);
0402     output.setY(vec.t()*_tensor[3][1]-vec.x()*_tensor[0][1]-
0403                 vec.y()*_tensor[1][1]-vec.z()*_tensor[2][1]);
0404     output.setZ(vec.t()*_tensor[3][2]-vec.x()*_tensor[0][2]-
0405                 vec.y()*_tensor[1][2]-vec.z()*_tensor[2][2]);
0406     output.setT(vec.t()*_tensor[3][3]-vec.x()*_tensor[0][3]-
0407                 vec.y()*_tensor[1][3]-vec.z()*_tensor[2][3]);
0408     return output;
0409   }
0410
0411   /**
0412    *  Second index dot product with polarization vector
0413    */
0414   template<typename ValueB>
0415   auto postDot(const LorentzVector<complex<ValueB> > & vec) const
0416     -> LorentzVector<decltype(vec.x()*this->xx())> {
0417     LorentzVector<decltype(vec.x()*this->xx())> output;
0418     output.setX(vec.t()*_tensor[0][3]-vec.x()*_tensor[0][0]-
0419                 vec.y()*_tensor[0][1]-vec.z()*_tensor[0][2]);
0420     output.setY(vec.t()*_tensor[1][3]-vec.x()*_tensor[1][0]-
0421                 vec.y()*_tensor[1][1]-vec.z()*_tensor[1][2]);
0422     output.setZ(vec.t()*_tensor[2][3]-vec.x()*_tensor[2][0]-
0423                 vec.y()*_tensor[2][1]-vec.z()*_tensor[2][2]);
0424     output.setT(vec.t()*_tensor[3][3]-vec.x()*_tensor[3][0]-
0425                 vec.y()*_tensor[3][1]-vec.z()*_tensor[3][2]);
0426     return output;
0427   }
0428
0429   /**
0430    *  First index dot product with momentum
0431    */
0432   auto preDot (const Lorentz5Momentum & vec) const
0433   -> LorentzVector<decltype(vec.x()*this->xx())>
0434   {
0435     LorentzVector<decltype(vec.x()*this->xx())> output;
0436     output.setX(vec.t()*_tensor[3][0]-vec.x()*_tensor[0][0]-
0437                 vec.y()*_tensor[1][0]-vec.z()*_tensor[2][0]);
0438     output.setY(vec.t()*_tensor[3][1]-vec.x()*_tensor[0][1]-
0439                 vec.y()*_tensor[1][1]-vec.z()*_tensor[2][1]);
0440     output.setZ(vec.t()*_tensor[3][2]-vec.x()*_tensor[0][2]-
0441                 vec.y()*_tensor[1][2]-vec.z()*_tensor[2][2]);
0442     output.setT(vec.t()*_tensor[3][3]-vec.x()*_tensor[0][3]-
0443                 vec.y()*_tensor[1][3]-vec.z()*_tensor[2][3]);
0444     return output;
0445   }
0446
0447   /**
0448    *  Second index dot product with momentum
0449    */
0450   auto postDot(const Lorentz5Momentum & vec) const
0451   -> LorentzVector<decltype(vec.x()*this->xx())>
0452   {
0453     LorentzVector<decltype(vec.x()*this->xx())> output;
0454     output.setX(vec.t()*_tensor[0][3]-vec.x()*_tensor[0][0]-
0455                 vec.y()*_tensor[0][1]-vec.z()*_tensor[0][2]);
0456     output.setY(vec.t()*_tensor[1][3]-vec.x()*_tensor[1][0]-
0457                 vec.y()*_tensor[1][1]-vec.z()*_tensor[1][2]);
0458     output.setZ(vec.t()*_tensor[2][3]-vec.x()*_tensor[2][0]-
0459                 vec.y()*_tensor[2][1]-vec.z()*_tensor[2][2]);
0460     output.setT(vec.t()*_tensor[3][3]-vec.x()*_tensor[3][0]-
0461                 vec.y()*_tensor[3][1]-vec.z()*_tensor[3][2]);
0462     return output;
0463   }
0464   //@}
0465 private:
0466
0467   /**
0468    * The components.
0469    */
0470   std::array<std::array<complex<Value>,4>,4> _tensor;
0471
0472 };
0473
0474 /**
0475  * Multiplication by a complex number.
0476  */
0477 template<typename T, typename U>
0478 inline auto
0479 operator*(complex<U> a, const LorentzTensor<T> & t)
0480 -> LorentzTensor<decltype(a.real()*t.xx().real())>
0481 {
0482   return
0483     {a*t.xx(), a*t.xy(), a*t.xz(), a*t.xt(),
0484      a*t.yx(), a*t.yy(), a*t.yz(), a*t.yt(),
0485      a*t.zx(), a*t.zy(), a*t.zz(), a*t.zt(),
0486      a*t.tx(), a*t.ty(), a*t.tz(), a*t.tt()};
0487 }
0488
0489 /**
0490  * Multiplication by a complex number.
0491  */
0492 template<typename T, typename U>
0493 inline auto
0494 operator*(const LorentzTensor<T> & t,complex<U> a)
0495 -> LorentzTensor<decltype(a.real()*t.xx().real())>
0496 {
0497   return
0498     {a*t.xx(), a*t.xy(), a*t.xz(), a*t.xt(),
0499      a*t.yx(), a*t.yy(), a*t.yz(), a*t.yt(),
0500      a*t.zx(), a*t.zy(), a*t.zz(), a*t.zt(),
0501      a*t.tx(), a*t.ty(), a*t.tz(), a*t.tt()};
0502 }
0503
0504 /**
0505  * Multiply a LorentzVector by a LorentzTensor.
0506  */
0507 template<typename T, typename U>
0508 inline auto
0509 operator*(const LorentzVector<U> & v,
0510           const LorentzTensor<T> & t)
0511 -> LorentzVector<decltype(v.t()*t(3,0))>
0512 {
0513   LorentzVector<decltype(v.t()*t(3,0))> outvec;
0514   outvec.setX( v.t()*t(3,0)-v.x()*t(0,0)
0515               -v.y()*t(1,0)-v.z()*t(2,0));
0516   outvec.setY( v.t()*t(3,1)-v.x()*t(0,1)
0517               -v.y()*t(1,1)-v.z()*t(2,1));
0518   outvec.setZ( v.t()*t(3,2)-v.x()*t(0,2)
0519               -v.y()*t(1,2)-v.z()*t(2,2));
0520   outvec.setT( v.t()*t(3,3)-v.x()*t(0,3)
0521               -v.y()*t(1,3)-v.z()*t(2,3));
0522   return outvec;
0523 }
0524
0525 /**
0526  * Multiply a LorentzTensor by a LorentzVector.
0527  */
0528 template<typename T, typename U>
0529 inline auto
0530 operator*(const LorentzTensor<T> & t, const LorentzVector<U> & v)
0531 -> LorentzVector<decltype(v.t()*t(0,3))>
0532 {
0533   LorentzVector<decltype(v.t()*t(0,3))> outvec;
0534   outvec.setX( v.t()*t(0,3)-v.x()*t(0,0)
0535               -v.y()*t(0,1)-v.z()*t(0,2));
0536   outvec.setY( v.t()*t(1,3)-v.x()*t(1,0)
0537               -v.y()*t(1,1)-v.z()*t(1,2));
0538   outvec.setZ( v.t()*t(2,3)-v.x()*t(2,0)
0539               -v.y()*t(2,1)-v.z()*t(2,2));
0540   outvec.setT( v.t()*t(3,3)-v.x()*t(3,0)
0541               -v.y()*t(3,1)-v.z()*t(3,2));
0542   return outvec;
0543 }
0544
0545 /**
0546  * Multiply a LorentzTensor by a LorentzTensor
0547  */
0548 template <typename T, typename U>
0549 inline auto
0550 operator*(const LorentzTensor<T> & t,
0551           const LorentzTensor<U> & u)
0552 -> decltype(t.xx()*u.xx())
0553 {
0554   using RetT = decltype(t.xx()*u.xx());
0555   RetT output=RetT(),temp;
0556   for(unsigned int ix=0;ix<4;++ix) {
0557     temp = t._tensor[ix][3]*u._tensor[ix][3];
0558     for(unsigned int iy=0;iy<3;++iy) {
0559       temp -= t._tensor[ix][iy]*u._tensor[ix][iy];
0560     }
0561     if(ix<3) output-=temp;
0562     else     output+=temp;
0563   }
0564   return output;
0565 }
0566
0567 }
0568 }
0569
0570 #ifndef ThePEG_TEMPLATES_IN_CC_FILE
0571 #include "LorentzTensor.tcc"
0572 #endif
0573
0574 #endif