ThePEG/Helicity/LorentzRank3Tensor.h

0001 // -*- C++ -*-
0002 //
0003 // LorentzRank3Tensor.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_LorentzRank3Tensor_H
0010 #define ThePEG_LorentzRank3Tensor_H
0011 // This is the declaration of the LorentzRank3Tensor class.
0012
0013 #include "ThePEG/Config/PhysicalQtyComplex.h"
0014 #include "ThePEG/Config/ThePEG.h"
0015 #include "LorentzTensor.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 LorentzRank3Tensor;
0022
0023 /**
0024  *  The LorentzRank3Tensor 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 LorentzRank3Tensor {
0038
0039 public:
0040
0041   /**
0042    * Default zero constructor.
0043    */
0044   LorentzRank3Tensor() = default;
0045
0046   /**
0047    * Get components by indices.
0048    */
0049   complex<Value> operator () (int i, int j, int k) const {
0050     assert( i>=0 && i<=3 && j>=0 && j<=3 && k>=0 && k<=3);
0051     return _tensor[i][j][k];
0052   }
0053
0054   /**
0055    * Set components by indices.
0056    */
0057   complex<Value> & operator () (int i, int j, int k) {
0058     assert( i>=0 && i<=3 && j>=0 && j<=3 && k>=0 && k<=3);
0059     return _tensor[i][j][k];
0060   }
0061   //@}
0062
0063   /** @name Transformations. */
0064   //@{
0065   /**
0066    * Standard Lorentz boost specifying the components of the beta vector.
0067    */
0068   LorentzRank3Tensor & boost(double,double,double);
0069
0070   /**
0071    * Standard Lorentz boost specifying the beta vector.
0072    */
0073   LorentzRank3Tensor<Value> & boost(const Boost & b) {
0074     return boost(b.x(), b.y(), b.z());
0075   }
0076
0077   /**
0078    * General Lorentz transformation
0079    */
0080   LorentzRank3Tensor & transform(const SpinOneLorentzRotation & r){
0081     unsigned int ix,iy,iz,ixa,iya,iza;
0082     LorentzRank3Tensor<Value> output;
0083     complex<Value> temp;
0084     for(ix=0;ix<4;++ix) {
0085       for(iy=0;iy<4;++iy) {
0086     for(iz=0;iz<4;++iz) {
0087       output(ix,iy,iz) = complex<Value>();
0088       for(ixa=0;ixa<4;++ixa) {
0089         for(iya=0;iya<4;++iya) {
0090           for(iza=0;iza<4;++iza)
0091          output(ix,iy,iz) += r(ix,ixa)*r(iy,iya)*r(iz,iza)*(*this)(ixa,iya,iza);
0092         }
0093       }
0094     }
0095       }
0096     }
0097     *this=output;
0098     return *this;
0099   }
0100
0101   /**
0102    * Return the complex conjugate.
0103    */
0104   LorentzRank3Tensor<Value> conjugate() {
0105     LorentzRank3Tensor<Value> output;
0106     for(unsigned int ix=0;ix<4;++ix) {
0107       for(unsigned int iy=0;iy<4;++iy) {
0108     for(unsigned int iz=0;iz<4;++iz) {
0109       output(ix,iy,iz) = conj(_tensor[ix][iy][iz]);
0110     }
0111       }
0112     }
0113     return output;
0114   }
0115
0116   //@}
0117
0118   /** @name Arithmetic operators. */
0119   //@{
0120   /**
0121    * Scaling with a complex number
0122    */
0123   LorentzRank3Tensor<Value> operator*=(Complex a) {
0124     for(int ix=0;ix<4;++ix)
0125       for(int iy=0;iy<4;++iy)
0126     for(int iz=0;iz<4;++iz) _tensor[ix][iy][iz]*=a;
0127     return *this;
0128   }
0129
0130   /**
0131    * Scalar product with other tensor
0132    */
0133   template <typename T, typename U>
0134   friend auto
0135   operator*(const LorentzRank3Tensor<T> & t, const LorentzRank3Tensor<U> & u) -> decltype(t.xx()*u.xx());
0136
0137   /**
0138    * Addition.
0139    */
0140   LorentzRank3Tensor<Value> operator+(const LorentzRank3Tensor<Value> & in) const {
0141     LorentzRank3Tensor<Value> output;
0142     for(int ix=0;ix<4;++ix)
0143       for(int iy=0;iy<4;++iy)
0144     for(int iz=0;iz<4;++iz) output(ix,iy,iz) = _tensor[ix][iy][iz] + in(ix,iy,iz);
0145   }
0146
0147   /**
0148    * Subtraction.
0149    */
0150   LorentzRank3Tensor<Value> operator-(const LorentzRank3Tensor<Value> & in) const {
0151     LorentzRank3Tensor<Value> output;
0152     for(int ix=0;ix<4;++ix)
0153       for(int iy=0;iy<4;++iy)
0154     for(int iz=0;iz<4;++iz) output(ix,iy,iz) = _tensor[ix][iy][iz] - in(ix,iy,iz);
0155   }
0156
0157   /**
0158    * Dot product with the ith index
0159    */
0160   template<typename ValueB>
0161   auto dot(const LorentzVector<complex<ValueB> > & vec, unsigned int iloc) const
0162     -> LorentzTensor<decltype(ValueB()*Value())> {
0163     LorentzTensor<decltype(ValueB()*Value())> output;
0164     if(iloc==0) {
0165       for(unsigned int iy=0;iy<4;++iy) {
0166     for(unsigned int iz=0;iz<4;++iz) {
0167       output(iy,iz) =
0168         vec.t()*_tensor[3][iy][iz] - vec.x()*_tensor[0][iy][iz] -
0169         vec.y()*_tensor[1][iy][iz] - vec.z()*_tensor[2][iy][iz];
0170     }
0171       }
0172     }
0173     else if(iloc==1) {
0174       for(unsigned int iy=0;iy<4;++iy) {
0175     for(unsigned int iz=0;iz<4;++iz) {
0176       output(iy,iz) =
0177         vec.t()*_tensor[iy][3][iz] - vec.x()*_tensor[iy][0][iz] -
0178         vec.y()*_tensor[iy][1][iz] - vec.z()*_tensor[iy][2][iz];
0179     }
0180       }
0181     }
0182     else if(iloc==2) {
0183       for(unsigned int iy=0;iy<4;++iy) {
0184     for(unsigned int iz=0;iz<4;++iz) {
0185       output(iy,iz) =
0186         vec.t()*_tensor[iy][iz][3] - vec.x()*_tensor[iy][iz][0] -
0187         vec.y()*_tensor[iy][iz][1] - vec.z()*_tensor[iy][iz][2];
0188     }
0189       }
0190     }
0191     else
0192       assert(false);
0193     return output;
0194   }
0195
0196   /**
0197    *  dot product with momentum
0198    */
0199   auto dot (const Lorentz5Momentum & vec,unsigned int iloc) const
0200   -> LorentzTensor<decltype(vec.x()*Value())>
0201   {
0202     LorentzTensor<decltype(vec.x()*Value())> output;
0203     if(iloc==0) {
0204       for(unsigned int iy=0;iy<4;++iy) {
0205     for(unsigned int iz=0;iz<4;++iz) {
0206       output(iy,iz) =
0207         vec.t()*_tensor[3][iy][iz] - vec.x()*_tensor[0][iy][iz] -
0208         vec.y()*_tensor[1][iy][iz] - vec.z()*_tensor[2][iy][iz];
0209     }
0210       }
0211     }
0212     else if(iloc==1) {
0213       for(unsigned int iy=0;iy<4;++iy) {
0214     for(unsigned int iz=0;iz<4;++iz) {
0215       output(iy,iz) =
0216         vec.t()*_tensor[iy][3][iz] - vec.x()*_tensor[iy][0][iz] -
0217         vec.y()*_tensor[iy][1][iz] - vec.z()*_tensor[iy][2][iz];
0218     }
0219       }
0220     }
0221     else if(iloc==2) {
0222       for(unsigned int iy=0;iy<4;++iy) {
0223     for(unsigned int iz=0;iz<4;++iz) {
0224       output(iy,iz) =
0225         vec.t()*_tensor[iy][iz][3] - vec.x()*_tensor[iy][iz][0] -
0226         vec.y()*_tensor[iy][iz][1] - vec.z()*_tensor[iy][iz][2];
0227     }
0228       }
0229     }
0230     else
0231       assert(false);
0232     return output;
0233   }
0234   //@}
0235
0236 private:
0237
0238   /**
0239    * The components.
0240    */
0241   std::array<std::array<std::array<complex<Value>,4>,4>,4> _tensor;
0242
0243 };
0244
0245 /**
0246  * Multiplication by a complex number.
0247  */
0248 template<typename T, typename U>
0249 inline auto
0250 operator*(complex<U> a, const LorentzRank3Tensor<T> & t) -> LorentzRank3Tensor<decltype(a.real()*t.xx().real())> {
0251   LorentzRank3Tensor<decltype(a.real()*t.xx().real())> output;
0252   for(int ix=0;ix<4;++ix)
0253     for(int iy=0;iy<4;++iy)
0254       for(int iz=0;iz<4;++iz) output(ix,iy,iz) = a*t(ix,iy,iz);
0255   return output;
0256 }
0257
0258 /**
0259  * Multiplication by a complex number.
0260  */
0261 template<typename T, typename U>
0262 inline auto
0263 operator*(const LorentzRank3Tensor<T> & t,complex<U> a) -> LorentzRank3Tensor<decltype(a.real()*t.xx().real())> {
0264   LorentzRank3Tensor<decltype(a.real()*t.xx().real())> output;
0265   for(int ix=0;ix<4;++ix)
0266     for(int iy=0;iy<4;++iy)
0267       for(int iz=0;iz<4;++iz) output(ix,iy,iz) = a*t(ix,iy,iz);
0268   return output;
0269 }
0270
0271 }
0272 }
0273
0274 #ifndef ThePEG_TEMPLATES_IN_CC_FILE
0275 #include "LorentzRank3Tensor.tcc"
0276 #endif
0277
0278 #endif