diff --git a/Helicity/LorentzTensor.cc b/Helicity/LorentzRank3Tensor.cc copy from Helicity/LorentzTensor.cc copy to Helicity/LorentzRank3Tensor.cc --- a/Helicity/LorentzTensor.cc +++ b/Helicity/LorentzRank3Tensor.cc @@ -1,12 +1,12 @@ // -*- C++ -*- // -// LorentzTensor.cc is a part of ThePEG - Toolkit for HEP Event Generation +// LorentzRank3Tensor.cc 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. // #ifdef ThePEG_TEMPLATES_IN_CC_FILE -#include "LorentzTensor.tcc" +#include "LorentzRank3Tensor.tcc" #include "ThePEG/Utilities/DescribeClass.h" #endif diff --git a/Helicity/LorentzTensor.fh b/Helicity/LorentzRank3Tensor.fh copy from Helicity/LorentzTensor.fh copy to Helicity/LorentzRank3Tensor.fh --- a/Helicity/LorentzTensor.fh +++ b/Helicity/LorentzRank3Tensor.fh @@ -1,17 +1,17 @@ // -*- C++ -*- // // This is the forward declaration of the FermionSpinInfo class. // -#ifndef ThePEG_LorentzTensor_FH -#define ThePEG_LorentzTensor_FH +#ifndef ThePEG_LorentzRank3Tensor_FH +#define ThePEG_LorentzRank3Tensor_FH namespace ThePEG { namespace Helicity { template -class LorentzTensor; +class LorentzRank3Tensor; } } #endif diff --git a/Helicity/LorentzTensor.h b/Helicity/LorentzRank3Tensor.h copy from Helicity/LorentzTensor.h copy to Helicity/LorentzRank3Tensor.h --- a/Helicity/LorentzTensor.h +++ b/Helicity/LorentzRank3Tensor.h @@ -1,533 +1,199 @@ // -*- C++ -*- // -// LorentzTensor.h is a part of ThePEG - Toolkit for HEP Event Generation +// 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_LorentzTensor_H -#define ThePEG_LorentzTensor_H -// This is the declaration of the LorentzTensor class. +#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 "LorentzPolarizationVector.h" namespace ThePEG { namespace Helicity { // compiler magic needs these pre-declarations to make friend templates work -template class LorentzTensor; +template class LorentzRank3Tensor; /** - * The LorentzTensor class is designed to implement the storage of a + * 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 -class LorentzTensor { +class LorentzRank3Tensor { public: - /** @name Standard constructors and destructors. */ - //@{ /** * Default zero constructor. */ - LorentzTensor() = default; - - /** - * Constructor specifyign all components. - */ - LorentzTensor(complex xx, complex xy, - complex xz, complex xt, - complex yx, complex yy, - complex yz, complex yt, - complex zx, complex zy, - complex zz, complex zt, - complex tx, complex ty, - complex tz, complex 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 xx() const {return _tensor[0][0];} - - /** - * Get y,x component. - */ - complex yx() const {return _tensor[1][0];} - /** - * Get z,x component. - */ - complex zx() const {return _tensor[2][0];} - - /** - * Get t,x component. - */ - complex tx() const {return _tensor[3][0];} - - /** - * Get x,y component. - */ - complex xy() const {return _tensor[0][1];} - - /** - * Get y,y component. - */ - complex yy() const {return _tensor[1][1];} - - /** - * Get z,y component. - */ - complex zy() const {return _tensor[2][1];} - - /** - * Get t,y component. - */ - complex ty() const {return _tensor[3][1];} - - /** - * Get x,z component. - */ - complex xz() const {return _tensor[0][2];} - - /** - * Get y,z component. - */ - complex yz() const {return _tensor[1][2];} - - /** - * Get z,z component. - */ - complex zz() const {return _tensor[2][2];} - - /** - * Get t,z component. - */ - complex tz() const {return _tensor[3][2];} - - /** - * Get x,t component. - */ - complex xt() const {return _tensor[0][3];} - - /** - * Get y,t component. - */ - complex yt() const {return _tensor[1][3];} - - /** - * Get z,t component. - */ - complex zt() const {return _tensor[2][3];} - - /** - * Get t,t component. - */ - complex tt() const {return _tensor[3][3];} - - /** - * Set x,x component. - */ - void setXX(complex a) {_tensor[0][0]=a;} - - /** - * Set y,x component. - */ - void setYX(complex a) {_tensor[1][0]=a;} - - /** - * Set z,x component. - */ - void setZX(complex a) {_tensor[2][0]=a;} - - /** - * Set t,x component. - */ - void setTX(complex a) {_tensor[3][0]=a;} - - /** - * Set x,y component. - */ - void setXY(complex a) {_tensor[0][1]=a;} - - /** - * Set y,y component. - */ - void setYY(complex a) {_tensor[1][1]=a;} - - /** - * Set z,y component. - */ - void setZY(complex a) {_tensor[2][1]=a;} - - /** - * Set t,y component. - */ - void setTY(complex a) {_tensor[3][1]=a;} - - /** - * Set x,z component. - */ - void setXZ(complex a) {_tensor[0][2]=a;} - - /** - * Set y,z component. - */ - void setYZ(complex a) {_tensor[1][2]=a;} - - /** - * Set z,z component. - */ - void setZZ(complex a) {_tensor[2][2]=a;} - - /** - * Set t,z component. - */ - void setTZ(complex a) {_tensor[3][2]=a;} - - /** - * Set x,t component. - */ - void setXT(complex a) {_tensor[0][3]=a;} - - /** - * Set y,t component. - */ - void setYT(complex a) {_tensor[1][3]=a;} - - /** - * Set z,t component. - */ - void setZT(complex a) {_tensor[2][3]=a;} - - /** - * Set t,t component. - */ - void setTT(complex a) {_tensor[3][3]=a;} + LorentzRank3Tensor() = default; /** * Get components by indices. */ - complex operator () (int i, int j) const { - assert( i>=0 && i<=3 && j>=0 && j<=3); - return _tensor[i][j]; + complex 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 & operator () (int i, int j) { - assert( i>=0 && i<=3 && j>=0 && j<=3); - return _tensor[i][j]; + complex & 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. */ - LorentzTensor & boost(double,double,double); + LorentzRank3Tensor & boost(double,double,double); /** * Standard Lorentz boost specifying the beta vector. */ - LorentzTensor & boost(const Boost & b) { + LorentzRank3Tensor & 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 output; + LorentzRank3Tensor & transform(const SpinOneLorentzRotation & r){ + unsigned int ix,iy,iz,ixa,iya,iza; + LorentzRank3Tensor output; complex temp; for(ix=0;ix<4;++ix) { for(iy=0;iy<4;++iy) { - temp=complex(); - 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; + for(iz=0;iz<4;++iz) { + output(ix,iy,iz) = complex(); + 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. */ - LorentzTensor conjugate() { - return LorentzTensor(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())); + LorentzRank3Tensor conjugate() { + LorentzRank3Tensor 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 */ - LorentzTensor operator*=(Complex a) { + LorentzRank3Tensor operator*=(Complex a) { for(int ix=0;ix<4;++ix) - for(int iy=0;iy<4;++iy) _tensor[ix][iy]*=a; + 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 friend auto - operator*(const LorentzTensor & t, const LorentzTensor & u) - -> decltype(t.xx()*u.xx()) - ; + operator*(const LorentzRank3Tensor & t, const LorentzRank3Tensor & u) -> decltype(t.xx()*u.xx()); /** * Addition. */ - LorentzTensor operator+(const LorentzTensor & in) const { - return LorentzTensor(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()); + LorentzRank3Tensor operator+(const LorentzRank3Tensor & in) const { + LorentzRank3Tensor 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. */ - LorentzTensor operator-(const LorentzTensor & in) const { - return LorentzTensor(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 trace() const { - return _tensor[3][3]-_tensor[0][0]-_tensor[1][1]-_tensor[2][2]; + LorentzRank3Tensor operator-(const LorentzRank3Tensor & in) const { + LorentzRank3Tensor 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); } //@} - /** - * Various dot products - */ - //@{ - /** - * First index dot product with polarization vector - */ - LorentzVector > preDot (const LorentzPolarizationVector & vec) const { - LorentzVector > 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 - */ - LorentzVector > postDot(const LorentzPolarizationVector & vec) const { - LorentzVector > 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 - -> LorentzVectorxx())> - { - LorentzVectorxx())> 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 - -> LorentzVectorxx())> - { - LorentzVectorxx())> 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,4>,4> _tensor; + std::array,4>,4>,4> _tensor; }; /** * Multiplication by a complex number. */ template inline auto -operator*(complex a, const LorentzTensor & t) --> LorentzTensor -{ - 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()}; +operator*(complex a, const LorentzRank3Tensor & t) -> LorentzRank3Tensor { + LorentzRank3Tensor 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 inline auto -operator*(const LorentzTensor & t,complex a) --> LorentzTensor -{ - 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 -inline auto -operator*(const LorentzVector & v, - const LorentzTensor & t) --> LorentzVector -{ - LorentzVector 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 -inline auto -operator*(const LorentzTensor & t, const LorentzVector & v) --> LorentzVector -{ - LorentzVector 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 -inline auto -operator*(const LorentzTensor & t, - const LorentzTensor & 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; - } +operator*(const LorentzRank3Tensor & t,complex a) -> LorentzRank3Tensor { + LorentzRank3Tensor 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 "LorentzTensor.tcc" +#include "LorentzRank3Tensor.tcc" #endif #endif diff --git a/Helicity/LorentzTensor.tcc b/Helicity/LorentzRank3Tensor.tcc copy from Helicity/LorentzTensor.tcc copy to Helicity/LorentzRank3Tensor.tcc --- a/Helicity/LorentzTensor.tcc +++ b/Helicity/LorentzRank3Tensor.tcc @@ -1,48 +1,53 @@ // -*- C++ -*- // -// LorentzTensor.tcc is a part of ThePEG - Toolkit for HEP Event Generation +// LorentzRank3Tensor.tcc 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. // namespace ThePEG { namespace Helicity { // general boost template -LorentzTensor & LorentzTensor::boost(double bx, double by, double bz) { +LorentzRank3Tensor & LorentzRank3Tensor::boost(double bx, double by, double bz) { // basic definitions double boostm[4][4]; double b2 = bx*bx+by*by+bz*bz; double gamma = 1.0/sqrt(1.0-b2); double gmmone = b2 >0 ? (gamma-1.)/b2 : 0.0; double vec[3]={bx,by,bz}; // compute the lorentz boost matrix for(unsigned int ix=0;ix<3;++ix) { for(unsigned int iy=0;iy<3;++iy){boostm[ix][iy]=vec[ix]*vec[iy]*gmmone;} boostm[ix][ix]+=1; boostm[ix][3]=gamma*vec[ix]; boostm[3][ix]=boostm[ix][3]; } boostm[3][3]=gamma; // apply the boost - LorentzTensor output; + LorentzRank3Tensor output; complex temp; - unsigned int ix,iy,ixa,iya; + unsigned int ix,iy,iz,ixa,iya,iza; for(ix=0;ix<4;++ix) { for(iy=0;iy<4;++iy) { - temp=0.; - for(ixa=0;ixa<4;++ixa) { - for(iya=0;iya<4;++iya) - {temp+=boostm[ix][ixa]*boostm[iy][iya]*(*this)(ixa,iya);} + for(iz=0;iz<4;++iz) { + temp=0.; + for(ixa=0;ixa<4;++ixa) { + for(iya=0;iya<4;++iya) { + for(iza=0;iza<4;++iza) { + temp+=boostm[ix][ixa]*boostm[iy][iya]*boostm[iz][iza]*(*this)(ixa,iya,iza); + } + } + } + output(ix,iy)=temp; } - output(ix,iy)=temp; } } *this=output; return *this; } } } diff --git a/Helicity/Makefile.am b/Helicity/Makefile.am --- a/Helicity/Makefile.am +++ b/Helicity/Makefile.am @@ -1,33 +1,34 @@ SUBDIRS = WaveFunction Vertex mySOURCES = LorentzSpinor.cc \ - LorentzSpinorBar.cc LorentzTensor.cc \ - FermionSpinInfo.cc VectorSpinInfo.cc ScalarSpinInfo.cc TensorSpinInfo.cc \ + LorentzSpinorBar.cc LorentzTensor.cc LorentzRank3Tensor.cc \ + FermionSpinInfo.cc VectorSpinInfo.cc ScalarSpinInfo.cc TensorSpinInfo.cc Rank3TensorSpinInfo.cc \ LorentzRSSpinor.cc LorentzRSSpinorBar.cc RSFermionSpinInfo.cc DOCFILES = LorentzPolarizationVector.h LorentzSpinor.h \ - LorentzSpinorBar.h LorentzTensor.h \ - FermionSpinInfo.h VectorSpinInfo.h ScalarSpinInfo.h TensorSpinInfo.h \ + LorentzSpinorBar.h LorentzTensor.h LorentzRank3Tensor.h \ + FermionSpinInfo.h VectorSpinInfo.h ScalarSpinInfo.h TensorSpinInfo.h Rank3TensorSpinInfo.h \ LorentzRSSpinor.h LorentzRSSpinorBar.h RSFermionSpinInfo.h \ epsilon.h HelicityFunctions.h INCLUDEFILES = $(DOCFILES) \ LorentzSpinor.fh LorentzSpinor.tcc LorentzSpinorBar.tcc \ LorentzSpinorBar.fh \ FermionSpinInfo.fh \ LorentzTensor.fh \ + LorentzRank3Tensor.fh \ VectorSpinInfo.fh \ ScalarSpinInfo.fh \ - TensorSpinInfo.fh LorentzTensor.tcc \ + TensorSpinInfo.fh Rank3TensorSpinInfo.fh \ + LorentzTensor.tcc LorentzRank3Tensor.tcc \ HelicityDefinitions.h \ LorentzRSSpinor.fh LorentzRSSpinorBar.fh RSFermionSpinInfo.fh \ LorentzRSSpinor.tcc LorentzRSSpinorBar.tcc noinst_LTLIBRARIES = libThePEGHelicity.la libThePEGHelicity_la_SOURCES = $(mySOURCES) $(INCLUDEFILES) libThePEGHelicity_la_LIBADD = Vertex/libThePEGVertex.la \ WaveFunction/libThePEGWaveFunction.la include $(top_srcdir)/Config/Makefile.aminclude - diff --git a/Helicity/TensorSpinInfo.cc b/Helicity/Rank3TensorSpinInfo.cc copy from Helicity/TensorSpinInfo.cc copy to Helicity/Rank3TensorSpinInfo.cc --- a/Helicity/TensorSpinInfo.cc +++ b/Helicity/Rank3TensorSpinInfo.cc @@ -1,41 +1,40 @@ // -*- C++ -*- // -// TensorSpinInfo.cc is a part of ThePEG - Toolkit for HEP Event Generation +// Rank3TensorSpinInfo.cc 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. // // // This is the implementation of the non-inlined, non-templated member -// functions of the TensorSpinInfo class. +// functions of the Rank3TensorSpinInfo class. // // Author: Peter Richardson // -#include "TensorSpinInfo.h" +#include "Rank3TensorSpinInfo.h" #include "ThePEG/Utilities/DescribeClass.h" using namespace ThePEG; using namespace ThePEG::Helicity; // The following static variable is needed for the type // description system in ThePEG. -DescribeNoPIOClass -describeThePEGTensorSpinInfo("ThePEG::TensorSpinInfo", "libThePEG.so"); +DescribeNoPIOClass +describeThePEGRank3TensorSpinInfo("ThePEG::Rank3TensorSpinInfo", "libThePEG.so"); -void TensorSpinInfo::Init() {} +void Rank3TensorSpinInfo::Init() {} -void TensorSpinInfo::transform(const LorentzMomentum & m, - const LorentzRotation & r) { +void Rank3TensorSpinInfo::transform(const LorentzMomentum & m, + const LorentzRotation & r) { if(isNear(m)) { - for(unsigned int ix=0;ix<5;++ix) _currentstates[ix].transform(r.one()); + for(unsigned int ix=0;ix<7;++ix) _currentstates[ix].transform(r.one()); SpinInfo::transform(m,r); } } -EIPtr TensorSpinInfo::clone() const -{ +EIPtr Rank3TensorSpinInfo::clone() const { tcSpinPtr temp=this; return const_ptr_cast(temp); } diff --git a/Helicity/TensorSpinInfo.fh b/Helicity/Rank3TensorSpinInfo.fh copy from Helicity/TensorSpinInfo.fh copy to Helicity/Rank3TensorSpinInfo.fh --- a/Helicity/TensorSpinInfo.fh +++ b/Helicity/Rank3TensorSpinInfo.fh @@ -1,18 +1,18 @@ // -*- C++ -*- // -// This is the forward declaration of the TensorSpinInfo class. +// This is the forward declaration of the Rank3TensorSpinInfo class. // -#ifndef ThePEG_TensorSpinInfo_FH -#define ThePEG_TensorSpinInfo_FH +#ifndef ThePEG_Rank3TensorSpinInfo_FH +#define ThePEG_Rank3TensorSpinInfo_FH #include "ThePEG/Config/Pointers.h" namespace ThePEG { namespace Helicity { -ThePEG_DECLARE_CLASS_POINTERS(TensorSpinInfo,TensorSpinPtr); +ThePEG_DECLARE_CLASS_POINTERS(Rank3TensorSpinInfo,Rank3TensorSpinPtr); } } #endif diff --git a/Helicity/TensorSpinInfo.h b/Helicity/Rank3TensorSpinInfo.h copy from Helicity/TensorSpinInfo.h copy to Helicity/Rank3TensorSpinInfo.h --- a/Helicity/TensorSpinInfo.h +++ b/Helicity/Rank3TensorSpinInfo.h @@ -1,194 +1,197 @@ // -*- C++ -*- // -// TensorSpinInfo.h is a part of ThePEG - Toolkit for HEP Event Generation +// Rank3TensorSpinInfo.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_TensorSpinInfo_H -#define THEPEG_TensorSpinInfo_H -// This is the declaration of the TensorSpinInfo class. +#ifndef THEPEG_Rank3TensorSpinInfo_H +#define THEPEG_Rank3TensorSpinInfo_H +// This is the declaration of the Rank3TensorSpinInfo class. #include "ThePEG/EventRecord/SpinInfo.h" -#include "ThePEG/Helicity/LorentzTensor.h" -#include "TensorSpinInfo.fh" -// #include "TensorSpinInfo.xh" +#include "ThePEG/Helicity/LorentzRank3Tensor.h" +#include "Rank3TensorSpinInfo.fh" +// #include "Rank3TensorSpinInfo.xh" #include namespace ThePEG { namespace Helicity { /** - * The TensorSpinInfo class is the implementation of the spin + * The Rank3TensorSpinInfo class is the implementation of the spin * information for tensor particles. It inherits from the SpinInfo * class and implements the storage of the basis tensors. * * These basis states should be set by either matrix elements or * decayers which are capable of generating spin correlation * information. * * The basis states in the rest frame of the particles can then be * accessed by decayers to produce the correct correlation. * - * N.B. in our convention 0 is the \f$-2\f$ helicity state, - * 1 is the \f$-1\f$ helicity state, - * 2 is the \f$0\f$ helicity state, - * 3 is the \f$+1\f$ helicity state and - * 4 is the \f$+2\f$ helicity state. + * N.B. in our convention + * 0 is the \f$-3\f$ helicity state, + * 1 is the \f$-2\f$ helicity state, + * 2 is the \f$-1\f$ helicity state, + * 3 is the \f$0\f$ helicity state, + * 4 is the \f$+1\f$ helicity state, + * 5 is the \f$+2\f$ helicity state and + * 6 is the \f$+3\f$ helicity state. * * @author Peter Richardson * */ -class TensorSpinInfo: public SpinInfo { +class Rank3TensorSpinInfo: public SpinInfo { public: /** @name Standard constructors and destructors. */ //@{ /** * Default constructor. */ - TensorSpinInfo() : SpinInfo(PDT::Spin2), _decaycalc(false) {} + Rank3TensorSpinInfo() : SpinInfo(PDT::Spin3), _decaycalc(false) {} /** * Standard Constructor. * @param p the production momentum. * @param time true if the particle is time-like. */ - TensorSpinInfo(const Lorentz5Momentum & p, bool time) - : SpinInfo(PDT::Spin2, p, time), _decaycalc(false) {} + Rank3TensorSpinInfo(const Lorentz5Momentum & p, bool time) + : SpinInfo(PDT::Spin3, p, time), _decaycalc(false) {} //@} public: /** @name Access the basis states. */ //@{ /** * Set the basis state, this is production state. * @param hel the helicity (0,1,2,3,4 as described above.) - * @param in the LorentzTensor for the given helicity. + * @param in the LorentzRank3Tensor for the given helicity. */ - void setBasisState(unsigned int hel, LorentzTensor in) const { - assert(hel<5); + void setBasisState(unsigned int hel, LorentzRank3Tensor in) const { + assert(hel<7); _productionstates[hel]=in; _currentstates [hel]=in; } /** * Set the basis state for the decay. * @param hel the helicity (0,1,2,3,4 as described above.) - * @param in the LorentzTensor for the given helicity. + * @param in the LorentzRank3Tensor for the given helicity. */ - void setDecayState(unsigned int hel, LorentzTensor in) const { - assert(hel<5); + void setDecayState(unsigned int hel, LorentzRank3Tensor in) const { + assert(hel<7); _decaycalc = true; _decaystates[hel] = in; } /** * Get the basis state for the production for the given helicity, \a * hel (0,1,2,3,4 as described above.) */ - const LorentzTensor & getProductionBasisState(unsigned int hel) const { - assert(hel<5); + const LorentzRank3Tensor & getProductionBasisState(unsigned int hel) const { + assert(hel<7); return _productionstates[hel]; } /** * Get the basis state for the current for the given helicity, \a * hel (0,1,2,3,4 as described above.) */ - const LorentzTensor & getCurrentBasisState(unsigned int hel) const { - assert(hel<5); + const LorentzRank3Tensor & getCurrentBasisState(unsigned int hel) const { + assert(hel<7); return _currentstates[hel]; } /** * Get the basis state for the decay for the given helicity, \a hel - * (0,1,2,3,4 as described above.) + * (0,1,2,3,4,5,6 as described above.) */ - const LorentzTensor & getDecayBasisState(unsigned int hel) const { - assert(hel<5); + const LorentzRank3Tensor & getDecayBasisState(unsigned int hel) const { + assert(hel<7); if(!_decaycalc) { - for(unsigned int ix=0;ix<5;++ix) + for(unsigned int ix=0;ix<7;++ix) _decaystates[ix]=_currentstates[ix].conjugate(); _decaycalc=true; } return _decaystates[hel]; } //@} /** * Perform a lorentz rotation of the spin information */ virtual void transform(const LorentzMomentum &,const LorentzRotation &); /** * Undecay */ virtual void undecay() const { _decaycalc=false; SpinInfo::undecay(); } /** * Reset */ virtual void reset() { undecay(); _currentstates = _productionstates; SpinInfo::reset(); } public: /** * Standard Init function. */ static void Init(); /** * Standard clone method. */ virtual EIPtr clone() const; private: /** * Private and non-existent assignment operator. */ - TensorSpinInfo & operator=(const TensorSpinInfo &) = delete; + Rank3TensorSpinInfo & operator=(const Rank3TensorSpinInfo &) = delete; private: /** * Basis states in the frame in which the particle was produced. */ - mutable std::array,5> _productionstates; + mutable std::array,7> _productionstates; /** * Basis states in the frame in which the particle decays. */ - mutable std::array,5> _decaystates; + mutable std::array,7> _decaystates; /** * Basis states in the current frame of the particle */ - mutable std::array,5> _currentstates; + mutable std::array,7> _currentstates; /** * True if the decay state has been set. */ mutable bool _decaycalc; }; } } namespace ThePEG { } -#endif /* THEPEG_TensorSpinInfo_H */ +#endif /* THEPEG_Rank3TensorSpinInfo_H */ diff --git a/Helicity/TensorSpinInfo.h b/Helicity/TensorSpinInfo.h --- a/Helicity/TensorSpinInfo.h +++ b/Helicity/TensorSpinInfo.h @@ -1,194 +1,190 @@ // -*- C++ -*- // // TensorSpinInfo.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_TensorSpinInfo_H #define THEPEG_TensorSpinInfo_H // This is the declaration of the TensorSpinInfo class. #include "ThePEG/EventRecord/SpinInfo.h" #include "ThePEG/Helicity/LorentzTensor.h" #include "TensorSpinInfo.fh" // #include "TensorSpinInfo.xh" #include namespace ThePEG { namespace Helicity { /** * The TensorSpinInfo class is the implementation of the spin * information for tensor particles. It inherits from the SpinInfo * class and implements the storage of the basis tensors. * * These basis states should be set by either matrix elements or * decayers which are capable of generating spin correlation * information. * * The basis states in the rest frame of the particles can then be * accessed by decayers to produce the correct correlation. * * N.B. in our convention 0 is the \f$-2\f$ helicity state, * 1 is the \f$-1\f$ helicity state, * 2 is the \f$0\f$ helicity state, * 3 is the \f$+1\f$ helicity state and * 4 is the \f$+2\f$ helicity state. * * @author Peter Richardson * */ class TensorSpinInfo: public SpinInfo { public: /** @name Standard constructors and destructors. */ //@{ /** * Default constructor. */ TensorSpinInfo() : SpinInfo(PDT::Spin2), _decaycalc(false) {} /** * Standard Constructor. * @param p the production momentum. * @param time true if the particle is time-like. */ TensorSpinInfo(const Lorentz5Momentum & p, bool time) : SpinInfo(PDT::Spin2, p, time), _decaycalc(false) {} //@} public: /** @name Access the basis states. */ //@{ /** * Set the basis state, this is production state. * @param hel the helicity (0,1,2,3,4 as described above.) * @param in the LorentzTensor for the given helicity. */ void setBasisState(unsigned int hel, LorentzTensor in) const { assert(hel<5); _productionstates[hel]=in; _currentstates [hel]=in; } /** * Set the basis state for the decay. * @param hel the helicity (0,1,2,3,4 as described above.) * @param in the LorentzTensor for the given helicity. */ void setDecayState(unsigned int hel, LorentzTensor in) const { assert(hel<5); _decaycalc = true; _decaystates[hel] = in; } /** * Get the basis state for the production for the given helicity, \a * hel (0,1,2,3,4 as described above.) */ const LorentzTensor & getProductionBasisState(unsigned int hel) const { assert(hel<5); return _productionstates[hel]; } /** * Get the basis state for the current for the given helicity, \a * hel (0,1,2,3,4 as described above.) */ const LorentzTensor & getCurrentBasisState(unsigned int hel) const { assert(hel<5); return _currentstates[hel]; } /** * Get the basis state for the decay for the given helicity, \a hel * (0,1,2,3,4 as described above.) */ const LorentzTensor & getDecayBasisState(unsigned int hel) const { assert(hel<5); if(!_decaycalc) { for(unsigned int ix=0;ix<5;++ix) _decaystates[ix]=_currentstates[ix].conjugate(); _decaycalc=true; } return _decaystates[hel]; } //@} /** * Perform a lorentz rotation of the spin information */ virtual void transform(const LorentzMomentum &,const LorentzRotation &); /** * Undecay */ virtual void undecay() const { _decaycalc=false; SpinInfo::undecay(); } /** * Reset */ virtual void reset() { undecay(); _currentstates = _productionstates; SpinInfo::reset(); } public: /** * Standard Init function. */ static void Init(); /** * Standard clone method. */ virtual EIPtr clone() const; private: /** * Private and non-existent assignment operator. */ TensorSpinInfo & operator=(const TensorSpinInfo &) = delete; private: /** * Basis states in the frame in which the particle was produced. */ mutable std::array,5> _productionstates; /** * Basis states in the frame in which the particle decays. */ mutable std::array,5> _decaystates; /** * Basis states in the current frame of the particle */ mutable std::array,5> _currentstates; /** * True if the decay state has been set. */ mutable bool _decaycalc; }; } } - -namespace ThePEG { - -} #endif /* THEPEG_TensorSpinInfo_H */ diff --git a/Helicity/WaveFunction/Makefile.am b/Helicity/WaveFunction/Makefile.am --- a/Helicity/WaveFunction/Makefile.am +++ b/Helicity/WaveFunction/Makefile.am @@ -1,10 +1,11 @@ noinst_LTLIBRARIES = libThePEGWaveFunction.la libThePEGWaveFunction_la_SOURCES = \ RSSpinorBarWaveFunction.cc RSSpinorBarWaveFunction.h \ RSSpinorWaveFunction.cc RSSpinorWaveFunction.h \ ScalarWaveFunction.h \ SpinorBarWaveFunction.cc SpinorBarWaveFunction.h\ SpinorWaveFunction.cc SpinorWaveFunction.h \ TensorWaveFunction.cc TensorWaveFunction.h \ +Rank3TensorWaveFunction.cc Rank3TensorWaveFunction.h \ VectorWaveFunction.cc VectorWaveFunction.h \ WaveFunctionBase.h diff --git a/Helicity/WaveFunction/TensorWaveFunction.cc b/Helicity/WaveFunction/Rank3TensorWaveFunction.cc copy from Helicity/WaveFunction/TensorWaveFunction.cc copy to Helicity/WaveFunction/Rank3TensorWaveFunction.cc --- a/Helicity/WaveFunction/TensorWaveFunction.cc +++ b/Helicity/WaveFunction/Rank3TensorWaveFunction.cc @@ -1,338 +1,356 @@ // -*- C++ -*- // -// TensorWaveFunction.cc is a part of ThePEG - Toolkit for HEP Event Generation +// Rank3TensorWaveFunction.cc 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. // // // This is the implementation of the non-inlined, non-templated member -// functions of the TensorWaveFunction class. +// functions of the Rank3TensorWaveFunction class. // // Author: Peter Richardson // -#include "TensorWaveFunction.h" +#include "Rank3TensorWaveFunction.h" using namespace ThePEG; using namespace ThePEG::Helicity; // calculate the actual wavefunction -void TensorWaveFunction::calculateWaveFunction(unsigned int ihel, TensorPhase tphase) { - int jhel=ihel-2; +void Rank3TensorWaveFunction::calculateWaveFunction(unsigned int ihel) { + int jhel=ihel-3; assert(direction()!=intermediate); // check for a valid helicty combination - assert((jhel<=2 && jhel>=-2 && mass() >ZERO) || - ((jhel==2 || jhel==-2) && mass()==ZERO)); + assert( (jhel<=3 && jhel>=-3 && mass() >ZERO) || + ((jhel==3 || jhel==-3) && mass()==ZERO)); // extract the momentum components double fact = direction()==outgoing ? -1. : 1; Energy ppx=fact*px(),ppy=fact*py(),ppz=fact*pz(),pee=fact*e(),pmm=mass(); // calculate some kinematic quantites; Energy2 pt2 = sqr(ppx)+sqr(ppy); Energy pabs = sqrt(pt2+sqr(ppz)); Energy pt = sqrt(pt2); // polarization vectors complex epsp[4],epsm[4],eps0[4]; // + helicity vector if needed - if(jhel>=0) { + if(jhel>=-1) { // calculate the overall phase - complexphase; - if(tphase==tensor_phase) { - phase = pt==ZERO ? 1. : complex(ppx/pt,-fact*ppy/pt); - } - else phase = 1.; - phase = phase*sqrt(0.5); + complex phase = pt==ZERO ? sqrt(0.5) : sqrt(0.5)*complex(ppx/pt,-fact*ppy/pt); // first the no pt case if(pt==ZERO) { double sgnz = ppz(0,-fact); epsp[2]=0.; epsp[3]=0.; } else { InvEnergy opabs=1./pabs; InvEnergy opt =1./pt; epsp[0]=phase*complex(-ppz*ppx*opabs*opt, fact*ppy*opt); epsp[1]=phase*complex(-ppz*ppy*opabs*opt, -fact*ppx*opt); epsp[2]=Complex(pt*opabs*phase); epsp[3]=0.; } } // - helicity vector if needed - if(jhel<=0) { + if(jhel<=1) { // calculate the overall phase - complex phase; - if(tphase==tensor_phase) { - phase = pt==ZERO ? 1. : complex(ppx/pt,fact*ppy/pt); - } - else phase = 1.; - phase = phase*sqrt(0.5); + complex phase = pt==ZERO ? sqrt(0.5) : sqrt(0.5)*complex(ppx/pt,fact*ppy/pt); // first the no pt case if(pt==ZERO) { - double sgnz; - if(ppz(0,-fact); epsm[2]=0.; epsm[3]=0.; } else { InvEnergy opabs=1./pabs; InvEnergy opt =1./pt; epsm[0]=phase*complex(ppz*ppx*opabs*opt, fact*ppy*opt); epsm[1]=phase*complex(ppz*ppy*opabs*opt, -fact*ppx*opt); epsm[2]=-Complex(pt*opabs*phase); epsm[3]=0.; } } // 0 helicity vector if needed - if(jhel<=1 && jhel>=-1) { + if(abs(jhel)!=3) { if(pabs==ZERO) { eps0[0] = 0.; eps0[1] = 0.; eps0[2] = 1.; eps0[3] = 0.; } else { InvEnergy empabs=pee/pmm/pabs; eps0[0] = empabs*ppx; eps0[1] = empabs*ppy; eps0[2] = empabs*ppz; eps0[3] = pabs/pmm; } } // put the polarization vectors together to get the wavefunction double ort; - switch (jhel) { + switch (jhel) { + case 3: + for(int ix=0;ix<4;++ix) + for(int iy=0;iy<4;++iy) + for(int iz=0;iz<4;++iz) _wf(ix,iy,iz)=epsp[ix]*epsp[iy]*epsp[iz]; + break; case 2: + ort = sqrt(0.5); for(int ix=0;ix<4;++ix) - for(int iy=0;iy<4;++iy) _wf(ix,iy)=epsp[ix]*epsp[iy]; + for(int iy=0;iy<4;++iy) + for(int iz=0;iz<4;++iz) + _wf(ix,iy,iz)= ort*(epsp[ix]*epsp[iy]*eps0[iz] + + epsp[ix]*eps0[iy]*epsp[iz] + + eps0[ix]*epsp[iy]*epsp[iz]); break; case 1: + ort = sqrt(1./15.); + for(int ix=0;ix<4;++ix) + for(int iy=0;iy<4;++iy) + for(int iz=0;iz<4;++iz) + _wf(ix,iy,iz)= ort*(epsp[ix]*epsp[iy]*epsm[iz] + + epsp[ix]*epsm[iy]*epsp[iz] + + epsm[ix]*epsp[iy]*epsp[iz] + +2.*(eps0[ix]*eps0[iy]*epsp[iz] + + eps0[ix]*epsp[iy]*eps0[iz] + + epsp[ix]*eps0[iy]*eps0[iz])); + break; + case 0: + ort = sqrt(0.1); + for(int ix=0;ix<4;++ix) + for(int iy=0;iy<4;++iy) + for(int iz=0;iz<4;++iz) + _wf(ix,iy,iz)= ort*(epsp[ix]*epsm[iy]*eps0[iz] + + epsp[ix]*eps0[iy]*epsm[iz] + + eps0[ix]*epsp[iy]*epsm[iz] + + eps0[ix]*epsm[iy]*epsp[iz] + + epsm[ix]*epsp[iy]*eps0[iz] + + epsm[ix]*eps0[iy]*epsp[iz] + + 2.*eps0[ix]*eps0[iy]*eps0[iz]); + break; + case -1: + ort = sqrt(1./15.); + for(int ix=0;ix<4;++ix) + for(int iy=0;iy<4;++iy) + for(int iz=0;iz<4;++iz) + _wf(ix,iy,iz)= ort*(epsm[ix]*epsm[iy]*epsp[iz] + + epsm[ix]*epsp[iy]*epsm[iz] + + epsp[ix]*epsm[iy]*epsm[iz] + +2.*(eps0[ix]*eps0[iy]*epsm[iz] + + eps0[ix]*epsm[iy]*eps0[iz] + + epsm[ix]*eps0[iy]*eps0[iz])); + break; + case -2: ort = sqrt(0.5); for(int ix=0;ix<4;++ix) - for(int iy=0;iy<4;++iy) _wf(ix,iy)=ort*( epsp[ix]*eps0[iy]+ - eps0[ix]*epsp[iy]); - break; - case 0: - ort = 1./sqrt(6.); + for(int iy=0;iy<4;++iy) + for(int iz=0;iz<4;++iz) + _wf(ix,iy,iz)= ort*(epsm[ix]*epsm[iy]*eps0[iz] + + epsm[ix]*eps0[iy]*epsm[iz] + + eps0[ix]*epsm[iy]*epsm[iz]); + break;; + case -3: for(int ix=0;ix<4;++ix) - for(int iy=0;iy<4;++iy) _wf(ix,iy)=ort*( epsp[ix]*epsm[iy] - + epsm[ix]*epsp[iy] - +2.*eps0[ix]*eps0[iy]); - break; - case -1: - ort = 1./sqrt(2.); - for(int ix=0;ix<4;++ix) - for(int iy=0;iy<4;++iy) _wf(ix,iy)=ort*( epsm[ix]*eps0[iy]+ - eps0[ix]*epsm[iy]); - break; - case -2: - for(int ix=0;ix<4;++ix) - for(int iy=0;iy<4;++iy) _wf(ix,iy)=epsm[ix]*epsm[iy]; + for(int iy=0;iy<4;++iy) + for(int iz=0;iz<4;++iz) _wf(ix,iy,iz)=epsm[ix]*epsm[iy]*epsm[iz]; break; default: ThePEG::Helicity::HelicityConsistencyError() - << "Invalid Helicity = " << jhel << " requested for Tensor" + << "Invalid Helicity = " << jhel << " requested for Rank3Tensor" << Exception::abortnow; break; } } -void TensorWaveFunction:: -calculateWaveFunctions(vector > & waves, - tPPtr particle,Direction dir, bool massless, - TensorPhase phase) { - tTensorSpinPtr inspin = !particle->spinInfo() ? tTensorSpinPtr() : - dynamic_ptr_cast(particle->spinInfo()); - waves.resize(5); +void Rank3TensorWaveFunction:: +calculateWaveFunctions(vector > & waves, + tPPtr particle,Direction dir, bool massless) { + tRank3TensorSpinPtr inspin = !particle->spinInfo() ? tRank3TensorSpinPtr() : + dynamic_ptr_cast(particle->spinInfo()); + waves.resize(7); if(inspin) { if(dir==outgoing) { - for(unsigned int ix=0;ix<5;++ix) + for(unsigned int ix=0;ix<7;++ix) waves[ix]=inspin->getProductionBasisState(ix); } else { inspin->decay(); - for(unsigned int ix=0;ix<5;++ix) + for(unsigned int ix=0;ix<7;++ix) waves[ix]=inspin->getDecayBasisState(ix); } } else { assert(!particle->spinInfo()); - TensorWaveFunction wave(particle->momentum(),particle->dataPtr(),0, - dir,phase); - for(unsigned int ix=0;ix<5;++ix) { - if(massless&&ix>0&&ix<5) { - waves[ix] = LorentzTensor(); + Rank3TensorWaveFunction wave(particle->momentum(),particle->dataPtr(),0,dir); + for(unsigned int ix=0;ix<7;++ix) { + if(massless&&ix>0&&ix<7) { + waves[ix] = LorentzRank3Tensor(); } else { if(ix!=0) wave.reset(ix); waves[ix] = wave.wave(); } } } } -void TensorWaveFunction:: -calculateWaveFunctions(vector & waves, - tPPtr particle, Direction dir, bool massless, - TensorPhase phase) { - tTensorSpinPtr inspin = !particle->spinInfo() ? tTensorSpinPtr() : - dynamic_ptr_cast(particle->spinInfo()); - waves.resize(5); +void Rank3TensorWaveFunction:: +calculateWaveFunctions(vector & waves, + tPPtr particle, Direction dir, bool massless) { + tRank3TensorSpinPtr inspin = !particle->spinInfo() ? tRank3TensorSpinPtr() : + dynamic_ptr_cast(particle->spinInfo()); + waves.resize(7); if(inspin) { if(dir==outgoing) { - for(unsigned int ix=0;ix<5;++ix) - waves[ix]=TensorWaveFunction(particle->momentum(), + for(unsigned int ix=0;ix<7;++ix) + waves[ix]=Rank3TensorWaveFunction(particle->momentum(), particle->dataPtr(), inspin->getProductionBasisState(ix),dir); } else { inspin->decay(); - for(unsigned int ix=0;ix<5;++ix) - waves[ix]=TensorWaveFunction(particle->momentum(), + for(unsigned int ix=0;ix<7;++ix) + waves[ix]=Rank3TensorWaveFunction(particle->momentum(), particle->dataPtr(), inspin->getDecayBasisState(ix),dir); } } else { assert(!particle->spinInfo()); - TensorWaveFunction wave(particle->momentum(),particle->dataPtr(),0, - dir,phase); - for(unsigned int ix=0;ix<5;++ix) { - if(massless&&ix>0&&ix<5) { - waves[ix] = TensorWaveFunction(particle->momentum(),particle->dataPtr(),dir); + Rank3TensorWaveFunction wave(particle->momentum(),particle->dataPtr(),0,dir); + for(unsigned int ix=0;ix<7;++ix) { + if(massless&&ix>0&&ix<7) { + waves[ix] = Rank3TensorWaveFunction(particle->momentum(),particle->dataPtr(),dir); } else { if(ix!=0) wave.reset(ix); waves[ix] = wave; } } } } -void TensorWaveFunction:: -calculateWaveFunctions(vector > & waves, +void Rank3TensorWaveFunction:: +calculateWaveFunctions(vector > & waves, RhoDMatrix & rho, - tPPtr particle,Direction dir,bool massless, - TensorPhase phase) { - tTensorSpinPtr inspin = !particle->spinInfo() ? tTensorSpinPtr() : - dynamic_ptr_cast(particle->spinInfo()); - waves.resize(5); + tPPtr particle,Direction dir,bool massless) { + tRank3TensorSpinPtr inspin = !particle->spinInfo() ? tRank3TensorSpinPtr() : + dynamic_ptr_cast(particle->spinInfo()); + waves.resize(7); if(inspin) { if(dir==outgoing) { - for(unsigned int ix=0;ix<5;++ix) + for(unsigned int ix=0;ix<7;++ix) waves[ix]=inspin->getProductionBasisState(ix); - rho = RhoDMatrix(PDT::Spin2); + rho = RhoDMatrix(PDT::Spin3); } else { inspin->decay(); - for(unsigned int ix=0;ix<5;++ix) + for(unsigned int ix=0;ix<7;++ix) waves[ix]=inspin->getDecayBasisState(ix); rho = inspin->rhoMatrix(); } } else { assert(!particle->spinInfo()); - TensorWaveFunction wave(particle->momentum(),particle->dataPtr(),0, - dir,phase); - for(unsigned int ix=0;ix<5;++ix) { - if(massless&&ix>0&&ix<5) { - waves[ix] = LorentzTensor(); + Rank3TensorWaveFunction wave(particle->momentum(),particle->dataPtr(),0,dir); + for(unsigned int ix=0;ix<7;++ix) { + if(massless&&ix>0&&ix<7) { + waves[ix] = LorentzRank3Tensor(); } else { if(ix!=0) wave.reset(ix); waves[ix] = wave.wave(); } } - rho = RhoDMatrix(PDT::Spin2); + rho = RhoDMatrix(PDT::Spin3); } } -void TensorWaveFunction:: -calculateWaveFunctions(vector & waves, +void Rank3TensorWaveFunction:: +calculateWaveFunctions(vector & waves, RhoDMatrix & rho, - tPPtr particle, Direction dir, bool massless, - TensorPhase phase) { - tTensorSpinPtr inspin = !particle->spinInfo() ? tTensorSpinPtr() : - dynamic_ptr_cast(particle->spinInfo()); - waves.resize(5); + tPPtr particle, Direction dir, bool massless) { + tRank3TensorSpinPtr inspin = !particle->spinInfo() ? tRank3TensorSpinPtr() : + dynamic_ptr_cast(particle->spinInfo()); + waves.resize(7); if(inspin) { if(dir==outgoing) { - for(unsigned int ix=0;ix<5;++ix) - waves[ix]=TensorWaveFunction(particle->momentum(), + for(unsigned int ix=0;ix<7;++ix) + waves[ix]=Rank3TensorWaveFunction(particle->momentum(), particle->dataPtr(), inspin->getProductionBasisState(ix),dir); - rho = RhoDMatrix(PDT::Spin2); + rho = RhoDMatrix(PDT::Spin3); } else { inspin->decay(); - for(unsigned int ix=0;ix<5;++ix) - waves[ix]=TensorWaveFunction(particle->momentum(), + for(unsigned int ix=0;ix<7;++ix) + waves[ix]=Rank3TensorWaveFunction(particle->momentum(), particle->dataPtr(), inspin->getDecayBasisState(ix),dir); rho = inspin->rhoMatrix(); } } else { assert(!particle->spinInfo()); - TensorWaveFunction wave(particle->momentum(),particle->dataPtr(),0, - dir,phase); - for(unsigned int ix=0;ix<5;++ix) { - if(massless&&ix>0&&ix<5) { - waves[ix] = TensorWaveFunction(particle->momentum(),particle->dataPtr(),dir); + Rank3TensorWaveFunction wave(particle->momentum(),particle->dataPtr(),0,dir); + for(unsigned int ix=0;ix<7;++ix) { + if(massless&&ix>0&&ix<7) { + waves[ix] = Rank3TensorWaveFunction(particle->momentum(),particle->dataPtr(),dir); } else { if(ix!=0) wave.reset(ix); waves[ix] = wave; } } - rho = RhoDMatrix(PDT::Spin2); + rho = RhoDMatrix(PDT::Spin3); } } -void TensorWaveFunction:: -constructSpinInfo(const vector > & waves, +void Rank3TensorWaveFunction:: +constructSpinInfo(const vector > & waves, tPPtr part,Direction dir, bool time,bool ) { - assert(waves.size()==5); - tTensorSpinPtr inspin = !part->spinInfo() ? tTensorSpinPtr() : - dynamic_ptr_cast(part->spinInfo()); + assert(waves.size()==7); + tRank3TensorSpinPtr inspin = !part->spinInfo() ? tRank3TensorSpinPtr() : + dynamic_ptr_cast(part->spinInfo()); if(inspin) { - for(unsigned int ix=0;ix<5;++ix) + for(unsigned int ix=0;ix<7;++ix) if(dir==outgoing) inspin->setBasisState(ix,waves[ix]); else inspin->setDecayState(ix,waves[ix]); } else { - TensorSpinPtr temp = new_ptr(TensorSpinInfo(part->momentum(),time)); + Rank3TensorSpinPtr temp = new_ptr(Rank3TensorSpinInfo(part->momentum(),time)); part->spinInfo(temp); - for(unsigned int ix=0;ix<5;++ix) + for(unsigned int ix=0;ix<7;++ix) if(dir==outgoing) temp->setBasisState(ix,waves[ix]); else temp->setDecayState(ix,waves[ix]); } } -void TensorWaveFunction:: -constructSpinInfo(const vector & waves, +void Rank3TensorWaveFunction:: +constructSpinInfo(const vector & waves, tPPtr part,Direction dir, bool time,bool ) { - assert(waves.size()==5); - tTensorSpinPtr inspin = !part->spinInfo() ? tTensorSpinPtr() : - dynamic_ptr_cast(part->spinInfo()); + assert(waves.size()==7); + tRank3TensorSpinPtr inspin = !part->spinInfo() ? tRank3TensorSpinPtr() : + dynamic_ptr_cast(part->spinInfo()); if(inspin) { - for(unsigned int ix=0;ix<5;++ix) + for(unsigned int ix=0;ix<7;++ix) if(dir==outgoing) inspin->setBasisState(ix,waves[ix].wave()); else inspin->setDecayState(ix,waves[ix].wave()); } else { - TensorSpinPtr temp = new_ptr(TensorSpinInfo(part->momentum(),time)); + Rank3TensorSpinPtr temp = new_ptr(Rank3TensorSpinInfo(part->momentum(),time)); part->spinInfo(temp); - for(unsigned int ix=0;ix<5;++ix) + for(unsigned int ix=0;ix<7;++ix) if(dir==outgoing) temp->setBasisState(ix,waves[ix].wave()); else temp->setDecayState(ix,waves[ix].wave()); } } diff --git a/Helicity/WaveFunction/TensorWaveFunction.h b/Helicity/WaveFunction/Rank3TensorWaveFunction.h copy from Helicity/WaveFunction/TensorWaveFunction.h copy to Helicity/WaveFunction/Rank3TensorWaveFunction.h --- a/Helicity/WaveFunction/TensorWaveFunction.h +++ b/Helicity/WaveFunction/Rank3TensorWaveFunction.h @@ -1,364 +1,223 @@ // -*- C++ -*- // -// TensorWaveFunction.h is a part of ThePEG - Toolkit for HEP Event Generation +// Rank3TensorWaveFunction.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_TensorWaveFunction_H -#define ThePEG_TensorWaveFunction_H +#ifndef ThePEG_Rank3TensorWaveFunction_H +#define ThePEG_Rank3TensorWaveFunction_H // -// This is the declaration of the TensorWaveFunction class. +// This is the declaration of the Rank3TensorWaveFunction class. // #include "WaveFunctionBase.h" #include "VectorWaveFunction.h" -#include -#include +#include +#include #include #include namespace ThePEG { namespace Helicity { -/**\ingroup Helicity - * Definition of the enumerated values of the phase to include in the - * calculation of the polarization tensor. - */ -enum TensorPhase { - tensor_phase, /**< Include the phase factor.*/ - tensor_nophase, /**< No phase-factor. */ - default_tensor_phase=tensor_nophase /**< Default option.*/ -}; - /** \ingroup Helicity * \author Peter Richardson * - * The TensorWaveFunction class is designed to store the wavefunction + * The Rank3TensorWaveFunction class is designed to store the wavefunction * of a tensor in a form suitable for use in helicity amplitude * calculations of the matrix element using a similar philosophy to the * FORTRAN HELAS code. * - * In addition to storing the tensor using the LorentzTensor class + * In addition to storing the tensor using the LorentzRank3Tensor class * it inherits from the WaveFunctionBase class to provide storage of * the momentum and ParticleData for the tensor particle. * * This class also contains the code which does the actually * calculation of the tensor wavefunction. * - * There are two choices available for the calculation of the - * wavefunction. These are set using the TensorPhase enumeration - * which specifies a default choice. - * The first choice, tensor_phase, includes a phase factor - * \f$\exp(\pm i \phi)\f$ for the \f$\pm\f$ helicity states while the second, - * tensor_nophase, does not. - * * N.B. In our convention - * 0 is the \f$-2\f$ helicity state, - * 1 is the \f$-1\f$ helicity state, - * 2 is the \f$ 0\f$ helicity state, - * 3 is the \f$+1\f$ helicity state and - * 4 is the \f$+2\f$ helicity state. + * 0 is the \f$-3\f$ helicity state, + * 1 is the \f$-2\f$ helicity state, + * 2 is the \f$-1\f$ helicity state, + * 3 is the \f$ 0\f$ helicity state, + * 4 is the \f$+1\f$ helicity state and + * 5 is the \f$+2\f$ helicity state. + * 6 is the \f$+3\f$ helicity state. * * @see WaveFunctionBase - * @see LorentzTensor + * @see LorentzRank3Tensor * @see VectorWaveFunction */ -class TensorWaveFunction : public WaveFunctionBase { +class Rank3TensorWaveFunction : public WaveFunctionBase { public: /** @name Standard constructors and destructors. */ //@{ /** * Constructor, set the momentum and Wavefunction, the direction can also * be specified. * @param p The momentum. * @param part The ParticleData pointer * @param wave The wavefunction, \e i.e. the polarization vector. * @param dir The direction of the particle. */ - TensorWaveFunction(const Lorentz5Momentum & p,tcPDPtr part, - const LorentzTensor & wave, - Direction dir=intermediate) - : WaveFunctionBase(p,part,dir), _wf(wave) - { - assert(iSpin()==PDT::Spin2); - } - - /** - * Constructor, set the momentum and the components of the tensor. - * @param p The momentum. - * @param part The ParticleData pointer - * @param xx The \f$xx\f$ component. - * @param xy The \f$xy\f$ component. - * @param xz The \f$xz\f$ component. - * @param xt The \f$xt\f$ component. - * @param yx The \f$yx\f$ component. - * @param yy The \f$yy\f$ component. - * @param yz The \f$yz\f$ component. - * @param yt The \f$yt\f$ component. - * @param zx The \f$zx\f$ component. - * @param zy The \f$zy\f$ component. - * @param zz The \f$zz\f$ component. - * @param zt The \f$zt\f$ component. - * @param tx The \f$tx\f$ component. - * @param ty The \f$ty\f$ component. - * @param tz The \f$tz\f$ component. - * @param tt The \f$tt\f$ component. - */ - TensorWaveFunction(const Lorentz5Momentum & p,tcPDPtr part, - Complex xx,Complex xy,Complex xz,Complex xt,Complex yx, - Complex yy,Complex yz,Complex yt,Complex zx,Complex zy, - Complex zz,Complex zt,Complex tx,Complex ty,Complex tz, - Complex tt) - : WaveFunctionBase(p,part), _wf(xx,xy,xz,xt, - yx,yy,yz,yt, - zx,zy,zz,zt, - tx,ty,tz,tt) - { - assert(iSpin()==PDT::Spin2); + Rank3TensorWaveFunction(const Lorentz5Momentum & p,tcPDPtr part, + const LorentzRank3Tensor & wave, + Direction dir=intermediate) + : WaveFunctionBase(p,part,dir), _wf(wave) { + assert(iSpin()==PDT::Spin3); } /** * Constructor, set the momentum, helicity, direction and optionally the phase * @param p The momentum. * @param part The ParticleData pointer * @param ihel The helicity (0,1,2,3,4 as described above.) * @param dir The direction. * @param phase The phase choice. */ - TensorWaveFunction(const Lorentz5Momentum & p,tcPDPtr part, - unsigned int ihel,Direction dir, - TensorPhase phase=default_tensor_phase) - : WaveFunctionBase(p,part,dir) - { - assert(iSpin()==PDT::Spin2); - calculateWaveFunction(ihel,phase); + Rank3TensorWaveFunction(const Lorentz5Momentum & p,tcPDPtr part, + unsigned int ihel,Direction dir) + : WaveFunctionBase(p,part,dir) { + assert(iSpin()==PDT::Spin3); + calculateWaveFunction(ihel); } /** * Constructor, set the 5-momentum and direction, zero the wavefunction. * @param p The momentum. * @param part The ParticleData pointer. * @param dir The direction. */ - TensorWaveFunction(const Lorentz5Momentum & p, - tcPDPtr part,Direction dir) - : WaveFunctionBase(p,part,dir), _wf() - { - assert(iSpin()==PDT::Spin2); + Rank3TensorWaveFunction(const Lorentz5Momentum & p, + tcPDPtr part,Direction dir) + : WaveFunctionBase(p,part,dir), _wf() { + assert(iSpin()==PDT::Spin3); } /** * Default constructor. */ - TensorWaveFunction() {} + Rank3TensorWaveFunction() {} /** * Special for spin correlations \todo make static? */ - TensorWaveFunction(vector & wave, - tPPtr part,Direction dir,bool time,bool massless, - bool=true, - TensorPhase phase=default_tensor_phase) { - calculateWaveFunctions(wave,part,dir,massless,phase); + Rank3TensorWaveFunction(vector & wave, + tPPtr part,Direction dir,bool time,bool massless, + bool=true) { + calculateWaveFunctions(wave,part,dir,massless); constructSpinInfo(wave,part,dir,time,massless); } //@} /** * Access to the wavefunction and its components. */ //@{ /** * Subscript operator for the wavefunction. */ - Complex operator ()(int i, int j) const { - return _wf(i,j); + Complex operator ()(int i, int j, int k) const { + return _wf(i,j,k); } /** * Set components by index. */ - Complex & operator () (int i, int j) { - return _wf(i,j); + Complex & operator () (int i, int j, int k) { + return _wf(i,j,k); } /** - * Return wavefunction as polarization vector. + * Return wavefunction as polarization rank-3 tensor . */ - const LorentzTensor & wave() const {return _wf;} - - /** - * Get the \f$xx\f$ component. - */ - Complex xx() const {return _wf.xx();} - - /** - * Get the \f$yx\f$ component. - */ - Complex yx() const {return _wf.yx();} - - /** - * Get the \f$zx\f$ component. - */ - Complex zx() const {return _wf.zx();} - - /** - * Get the \f$tx\f$ component. - */ - Complex tx() const {return _wf.tx();} - - /** - * Get the \f$xy\f$ component. - */ - Complex xy() const {return _wf.xy();} - - /** - * Get the \f$yy\f$ component. - */ - Complex yy() const {return _wf.yy();} - - /** - * Get the \f$zy\f$ component. - */ - Complex zy() const {return _wf.zy();} - - /** - * Get the \f$ty\f$ component. - */ - Complex ty() const {return _wf.ty();} - - /** - * Get the \f$xz\f$ component. - */ - Complex xz() const {return _wf.xz();} - - /** - * Get the \f$yz\f$ component. - */ - Complex yz() const {return _wf.yz();} - - /** - * Get the \f$zz\f$ component. - */ - Complex zz() const {return _wf.zz();} - - /** - * Get the \f$tz\f$ component. - */ - Complex tz() const {return _wf.tz();} - - /** - * Get the \f$xt\f$ component. - */ - Complex xt() const {return _wf.xt();} - - /** - * Get the \f$yt\f$ component. - */ - Complex yt() const {return _wf.yt();} - - /** - * Get the \f$zt\f$ component. - */ - Complex zt() const {return _wf.zt();} - - /** - * Get the \f$tt\f$ component. - */ - Complex tt() const {return _wf.tt();} + const LorentzRank3Tensor & wave() const {return _wf;} //@} /** * Reset functions. */ //@{ /** * Reset helicity (recalculate the tensor ). * @param ihel The new helicity (0,1,2,3,4 as described above.) - * @param phase The phase choice. */ - void reset(unsigned int ihel,TensorPhase phase=default_tensor_phase) { - calculateWaveFunction(ihel,phase); + void reset(unsigned int ihel) { + calculateWaveFunction(ihel); } //@} public: /** * Perform the Lorentz transformation of the wave function */ void transform(const LorentzRotation & r) { _wf.transform(r); transformMomentum(r); } public: /** * Calculate the wavefunctions */ - static void calculateWaveFunctions(vector > & waves, - tPPtr particle,Direction,bool massless, - TensorPhase phase=default_tensor_phase); + static void calculateWaveFunctions(vector > & waves, + tPPtr particle,Direction,bool massless); /** * Calculate the wavefunctions */ - static void calculateWaveFunctions(vector & waves, - tPPtr particle,Direction,bool massless, - TensorPhase phase=default_tensor_phase); + static void calculateWaveFunctions(vector & waves, + tPPtr particle,Direction,bool massless); /** * Calculate the wavefunctions */ - static void calculateWaveFunctions(vector > & waves, + static void calculateWaveFunctions(vector > & waves, RhoDMatrix & rho, - tPPtr particle,Direction,bool massless, - TensorPhase phase=default_tensor_phase); + tPPtr particle,Direction,bool massless); /** * Calculate the wavefunctions */ - static void calculateWaveFunctions(vector & waves, + static void calculateWaveFunctions(vector & waves, RhoDMatrix & rho, - tPPtr particle,Direction,bool massless, - TensorPhase phase=default_tensor_phase); + tPPtr particle,Direction,bool massless); /** * Construct the SpinInfo object */ - static void constructSpinInfo(const vector > & waves, + static void constructSpinInfo(const vector > & waves, tPPtr part,Direction dir, bool time,bool massless); /** * Construct the SpinInfo object */ - static void constructSpinInfo(const vector & waves, + static void constructSpinInfo(const vector & waves, tPPtr part,Direction dir, bool time,bool massless); private: /** * Calculate the wavefunction. * @param ihel The helicity (0,1,2,3,4 as described above.) - * @param phase The phase choice. */ - void calculateWaveFunction(unsigned int ihel, - TensorPhase phase=default_tensor_phase); + void calculateWaveFunction(unsigned int ihel); private: /** - * Storage of the wavefunction as a Lorentz Tensor. + * Storage of the wavefunction as a Lorentz Rank3Tensor. */ - LorentzTensor _wf; + LorentzRank3Tensor _wf; }; } } #endif diff --git a/Vectors/SpinOneLorentzRotation.h b/Vectors/SpinOneLorentzRotation.h --- a/Vectors/SpinOneLorentzRotation.h +++ b/Vectors/SpinOneLorentzRotation.h @@ -1,394 +1,396 @@ // -*- C++ -*- // // SpinOneLorentzRotation.h is a part of ThePEG - Toolkit for HEP Event Generation // Copyright (C) 1999-2019 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_SpinOneLorentzRotation_H #define ThePEG_SpinOneLorentzRotation_H #include "ThePEG/Helicity/HelicityDefinitions.h" +#include "ThePEG/Helicity/LorentzRank3Tensor.fh" #include "ThePEG/Helicity/LorentzTensor.fh" #include "ThePEG/Helicity/LorentzRSSpinor.fh" #include "ThePEG/Helicity/LorentzRSSpinorBar.fh" #include "ThreeVector.h" #include namespace ThePEG { /** * The SpinOneLorentzRotation class is ... */ class SpinOneLorentzRotation { public: /** @name Constructors and destructor. */ //@{ /** * Default constructor. Gives a unit matrix. */ SpinOneLorentzRotation() { xx_() = yy_() = zz_() = tt_() = 1.0; } /** * Constructor giving the components of a Lorentz boost. * @param bx The x-component of the boost * @param by The y-component of the boost * @param bz The z-component of the boost * @param gamma The \f$\gamma\f$ factor (optional) */ SpinOneLorentzRotation (double bx, double by, double bz, double gamma=-1.) { setBoost(bx,by,bz,gamma); } /** * Constructor giving the vector for a Lorentz boost. * @param b The boost vector * @param gamma The \f$\gamma\f$ factor (optional) */ explicit SpinOneLorentzRotation (const Boost & b, double gamma=-1.) { setBoost(b.x(), b.y(), b.z(),gamma); } //@} /** * Returns true if the Identity matrix. */ bool isIdentity() const; /** * Return the inverse. */ SpinOneLorentzRotation inverse() const; /** * Inverts the SpinOneLorentzRotation matrix. */ SpinOneLorentzRotation & invert() { return *this = inverse(); } /** * output operator */ std::ostream & print( std::ostream & os ) const; /** @name Set methods for speical cases of simple rotations and boosts */ //@{ /** * Specify the components of a Lorentz Boost * @param bx The x-component of the boost * @param by The y-component of the boost * @param bz The z-component of the boost * @param gamma The \f$\gamma\f$ factor (optional) */ SpinOneLorentzRotation & setBoost (double bx, double by, double bz, double gamma=-1.); /** * Specify a Lorentz Boost as a vector * @param b The boost vector * @param gamma The \f$\gamma\f$ factor (optional) */ SpinOneLorentzRotation & setBoost (const Boost & b, double gamma=-1.) { return setBoost(b.x(), b.y(), b.z(),gamma); } /** * Specify a rotation about a general axis by the angle given. * @param delta The angle * @param axis The axis */ SpinOneLorentzRotation & setRotate(double delta, const Axis & axis); /** * Specify a rotation by the given angle about the x-axis * @param angle The rotation angle */ SpinOneLorentzRotation & setRotateX (double angle); /** * Specify a rotation by the given angle about the y-axis * @param angle The rotation angle */ SpinOneLorentzRotation & setRotateY (double angle); /** * Specify a rotation by the given angle about the z-axis * @param angle The rotation angle */ SpinOneLorentzRotation & setRotateZ (double angle); //@} /** @name Access methods for the components of the spin-1 rotation */ //@{ /** * The xx component */ double xx() const { return matrix_[ 0]; } /** * The xy component */ double xy() const { return matrix_[ 1]; } /** * The xz component */ double xz() const { return matrix_[ 2]; } /** * The xt component */ double xt() const { return matrix_[ 3]; } /** * The yx component */ double yx() const { return matrix_[ 4]; } /** * The yy component */ double yy() const { return matrix_[ 5]; } /** * The yz component */ double yz() const { return matrix_[ 6]; } /** * The yt component */ double yt() const { return matrix_[ 7]; } /** * The zx component */ double zx() const { return matrix_[ 8]; } /** * The zy component */ double zy() const { return matrix_[ 9]; } /** * The zz component */ double zz() const { return matrix_[10]; } /** * The zt component */ double zt() const { return matrix_[11]; } /** * The tx component */ double tx() const { return matrix_[12]; } /** * The ty component */ double ty() const { return matrix_[13]; } /** * The tz component */ double tz() const { return matrix_[14]; } /** * The tt component */ double tt() const { return matrix_[15]; } //@} /** @name Transformation and product members */ //@{ /** * Product with a LorentzVector simply returns the rotated vector. */ template LorentzVector operator*(const LorentzVector & v) const { return LorentzVector (xx()*v.x() + xy()*v.y() + xz()*v.z() + xt()*v.t(), yx()*v.x() + yy()*v.y() + yz()*v.z() + yt()*v.t(), zx()*v.x() + zy()*v.y() + zz()*v.z() + zt()*v.t(), tx()*v.x() + ty()*v.y() + tz()*v.z() + tt()*v.t()); } /** * Product with a Lorentz5Vector simply returns the rotated vector. */ template Lorentz5Vector operator*(const Lorentz5Vector & v) const { return Lorentz5Vector (xx()*v.x() + xy()*v.y() + xz()*v.z() + xt()*v.t(), yx()*v.x() + yy()*v.y() + yz()*v.z() + yt()*v.t(), zx()*v.x() + zy()*v.y() + zz()*v.z() + zt()*v.t(), tx()*v.x() + ty()*v.y() + tz()*v.z() + tt()*v.t()); } /** * Product of two LorentzRotations (this) * lt - matrix multiplication * @param lt The LorentzRotation we are multiplying */ SpinOneLorentzRotation operator * (const SpinOneLorentzRotation & lt) const; /** * Multiply by and assign a*=b becomes a= a*b */ SpinOneLorentzRotation & operator *= (const SpinOneLorentzRotation & lt) { return *this = *this * lt; } /** * Transform (similar to *= but a.transform(b) becomes a = b*a */ SpinOneLorentzRotation & transform (const SpinOneLorentzRotation & lt) { return *this = lt * *this; } /** * Rotation around the x-axis; equivalent to LT = RotationX(delta) * LT */ SpinOneLorentzRotation & rotateX(double delta) { SpinOneLorentzRotation tmp; tmp.setRotateX(delta); return *this = tmp * *this; } /** * Rotation around the y-axis; equivalent to LT = RotationY(delta) * LT */ SpinOneLorentzRotation & rotateY(double delta) { SpinOneLorentzRotation tmp; tmp.setRotateY(delta); return *this = tmp * *this; } /** * Rotation around the z-axis; equivalent to LT = RotationZ(delta) * LT */ SpinOneLorentzRotation & rotateZ(double delta) { SpinOneLorentzRotation tmp; tmp.setRotateZ(delta); return *this = tmp * *this; } /** * Rotation around specified vector - LT = Rotation(delta,axis)*LT */ SpinOneLorentzRotation & rotate(double delta, const Axis & axis) { SpinOneLorentzRotation tmp; tmp.setRotate(delta, axis); return *this = tmp * *this; } /** * Pure boost along the x-axis; equivalent to LT = BoostX(beta) * LT */ SpinOneLorentzRotation & boostX(double beta) { return *this = SpinOneLorentzRotation(beta,0,0) * *this; } /** * Pure boost along the y-axis; equivalent to LT = BoostX(beta) * LT */ SpinOneLorentzRotation & boostY(double beta) { return *this = SpinOneLorentzRotation(0,beta,0) * *this; } /** * Pure boost along the z-axis; equivalent to LT = BoostX(beta) * LT */ SpinOneLorentzRotation & boostZ(double beta) { return *this = SpinOneLorentzRotation(0,0,beta) * *this; } /** * boost equivalent to LT = Boost(bx,by,bz) * LT * @param bx The x-component of the boost * @param by The y-component of the boost * @param bz The z-component of the boost * @param gamma The \f$\gamma\f$ factor (optional) */ SpinOneLorentzRotation & boost(double bx, double by, double bz, double gamma=-1.) { return *this = SpinOneLorentzRotation(bx,by,bz,gamma) * *this; } /** * boost equivalent to LT = Boost(bv) * LT * @param b The boost vector * @param gamma The \f$\gamma\f$ factor (optional) */ SpinOneLorentzRotation & boost(const Boost & b, double gamma=-1.) { return *this = SpinOneLorentzRotation(b.x(),b.y(),b.z(),gamma) * *this; } //@} private: template friend class Helicity::LorentzTensor; + template friend class Helicity::LorentzRank3Tensor; template friend class Helicity::LorentzRSSpinor; template friend class Helicity::LorentzRSSpinorBar; /// Matrix components, order: \f$(xx, xy, \ldots, tz, tt)\f$. array matrix_ = {}; /// Constructor from doubles. SpinOneLorentzRotation (double xx, double xy, double xz, double xt, double yx, double yy, double yz, double yt, double zx, double zy, double zz, double zt, double tx, double ty, double tz, double tt); /// Component access by index: x=0, t=3. double operator()(unsigned int i, unsigned int j) const { return matrix_[4*i + j]; } /// @name Component access. //@{ double & xx_() { return matrix_[ 0]; } double & xy_() { return matrix_[ 1]; } double & xz_() { return matrix_[ 2]; } double & xt_() { return matrix_[ 3]; } double & yx_() { return matrix_[ 4]; } double & yy_() { return matrix_[ 5]; } double & yz_() { return matrix_[ 6]; } double & yt_() { return matrix_[ 7]; } double & zx_() { return matrix_[ 8]; } double & zy_() { return matrix_[ 9]; } double & zz_() { return matrix_[10]; } double & zt_() { return matrix_[11]; } double & tx_() { return matrix_[12]; } double & ty_() { return matrix_[13]; } double & tz_() { return matrix_[14]; } double & tt_() { return matrix_[15]; } //@} }; /** * output operator */ inline std::ostream & operator<< ( std::ostream & os, const SpinOneLorentzRotation& lt ) { return lt.print(os); } } #endif /* ThePEG_SpinOneLorentzRotation_H */