#ifndef ATOOLS_Math_Vec4_H #define ATOOLS_Math_Vec4_H #include "ATOOLS/Math/MathTools.H" namespace ATOOLS { template class Vec4; template class Vec3; template class Vec4 { Scalar m_x[4]; template friend class Vec4; public: inline Vec4() { m_x[0]=m_x[1]=m_x[2]=m_x[3]=Scalar(0.0); } template inline Vec4(const Vec4 & vec) { m_x[0] = vec[0]; m_x[1] = vec[1]; m_x[2] = vec[2]; m_x[3] = vec[3]; } inline Vec4(const Scalar& x0, const Scalar& x1, const Scalar& x2, const Scalar& x3) { m_x[0] = x0; m_x[1] = x1; m_x[2] = x2; m_x[3] = x3; } inline Vec4(const Scalar& E, const Vec3& xn) { m_x[0] = E; m_x[1] = xn[1]; m_x[2] = xn[2]; m_x[3] = xn[3]; } inline Scalar& operator[] (int i) { return m_x[i]; } inline const Scalar operator[] (int i) const { return m_x[i]; } inline Vec4& operator+=(const Vec4& v) { m_x[0] += v[0]; m_x[1] += v[1]; m_x[2] += v[2]; m_x[3] += v[3]; return *this; } inline Vec4& operator-=(const Vec4& v) { m_x[0] -= v[0]; m_x[1] -= v[1]; m_x[2] -= v[2]; m_x[3] -= v[3]; return *this; } inline Vec4& operator*= (const Scalar& scal) { m_x[0] *= scal; m_x[1] *= scal; m_x[2] *= scal; m_x[3] *= scal; return *this; } inline Vec4 operator-() const { return Vec4(-m_x[0],-m_x[1],-m_x[2],-m_x[3]); } template inline PROMOTE(Scalar,Scalar2) operator*(const Vec4& v2) const { return m_x[0]*v2[0]-m_x[1]*v2[1]-m_x[2]*v2[2]-m_x[3]*v2[3]; } template inline Vec4 operator+ (const Vec4& v2) const { return Vec4 (m_x[0]+v2[0],m_x[1]+v2[1],m_x[2]+v2[2],m_x[3]+v2[3]); } template inline Vec4 operator- (const Vec4& v2) const { return Vec4 (m_x[0]-v2[0],m_x[1]-v2[1],m_x[2]-v2[2],m_x[3]-v2[3]); } template inline Vec4 operator* (const Scalar2& s) const { return Vec4 (s*m_x[0],s*m_x[1],s*m_x[2],s*m_x[3]); } template inline Vec4 operator/ (const Scalar2& scal) const { return (*this)*(1./scal); } inline bool Nan() const { for(short unsigned int i(0);i<4;++i) if (ATOOLS::IsNan(m_x[i])) return true; return false; } inline bool IsZero() const { for(short unsigned int i(0);i<4;++i) if (!ATOOLS::IsZero(m_x[i])) return false; return true; } Scalar Abs2() const { return (m_x[0]+m_x[3])*(m_x[0]-m_x[3])-m_x[1]*m_x[1]-m_x[2]*m_x[2]; } Scalar RelAbs2() const { const auto abs2 = Abs2(); return (abs2 == 0) ? 0 : abs2/(m_x[0]*m_x[0]); } Scalar Abs() const { return sqrt(Abs2()); } double Mass() const { return sqrt(ATOOLS::Abs(Abs2())); } Scalar P() const { return PSpat(); } Scalar PPerp2() const { return m_x[1]*m_x[1]+m_x[2]*m_x[2]; } Scalar PPerp() const { return sqrt(PPerp2()); } Scalar MPerp2() const { return (m_x[0]+m_x[3])*(m_x[0]-m_x[3]); } Scalar MPerp() const { return sqrt(MPerp2()); } Scalar MPerp2(const Vec4 &ref) const { return Abs2()+PPerp2(ref); } Scalar MPerp(const Vec4 &ref) const { return sqrt(MPerp2(ref)); } Scalar EPerp2() const { return m_x[0]*m_x[0]*PPerp2()/PSpat2(); } Scalar EPerp() const { return sqrt(EPerp2()); } Scalar PPlus() const { return m_x[0]+m_x[3]; } Scalar PMinus() const { return m_x[0]-m_x[3]; } Scalar PSpat2() const { return m_x[1]*m_x[1]+m_x[2]*m_x[2]+m_x[3]*m_x[3]; } Scalar PSpat() const { return sqrt(PSpat2()); } Scalar Y() const { return 0.5*log(PPlus()/PMinus()); } inline Vec4 Perp() const { return Vec4(0.,m_x[1],m_x[2],0.); } inline Vec4 Long() const { return Vec4(m_x[0],0.,0.,m_x[3]); } inline Vec4 Plus() const { Scalar pplus=0.5*PPlus(); return Vec4(pplus,0.0,0.0,pplus); } inline Vec4 Minus() const { Scalar pminus=0.5*PMinus(); return Vec4(pminus,0.0,0.0,-pminus); } Scalar PPerp2(const Vec4& ref) const { Scalar ppref = PPerp(ref); return ppref*ppref; } Scalar PPerp(const Vec4& ref) const { Vec3 perp=1./ATOOLS::Max(Vec3(ref).Abs(),1.e-12)* Vec3(ref); perp=Vec3(*this)-perp*(perp*Vec3(*this)); return perp.Abs(); } // some functions which only make sense for doubles // and thus have to be specialised for the type Scalar CosPhi() const; Scalar SinPhi() const; Scalar Phi() const; Scalar CosTheta() const; Scalar SinTheta() const; Scalar Theta() const; Scalar Eta() const; Scalar CosTheta(const Vec4& ref) const; Scalar Theta(const Vec4& ref) const; Scalar Eta(const Vec4& ref) const; Scalar CosDPhi(const Vec4& ref) const; Scalar DPhi(const Vec4& ref) const; Scalar DEta(const Vec4& ref) const; Scalar DY(const Vec4& ref) const; Scalar DR(const Vec4& ref) const; Scalar DR2(const Vec4& ref) const; Scalar DRy(const Vec4& ref) const; Scalar DR2y(const Vec4& ref) const; // standard vectors: const static Vec4 XVEC; const static Vec4 YVEC; const static Vec4 ZVEC; }; template inline Vec4 operator* (const Scal2& s, const Vec4& v) { return Vec4(v*s); } template inline Vec4 cross(const Vec4&q, const Vec4&r, const Vec4&s) { PROMOTE(S1,PROMOTE(S2,S3)) r0s1mr1s0(r[0]*s[1]-r[1]*s[0]), r0s2mr2s0(r[0]*s[2]-r[2]*s[0]), r0s3mr3s0(r[0]*s[3]-r[3]*s[0]), r1s2mr2s1(r[1]*s[2]-r[2]*s[1]), r1s3mr3s1(r[1]*s[3]-r[3]*s[1]), r2s3mr3s2(r[2]*s[3]-r[3]*s[2]); return Vec4 (-q[1]*r2s3mr3s2+q[2]*r1s3mr3s1-q[3]*r1s2mr2s1, -q[0]*r2s3mr3s2+q[2]*r0s3mr3s0-q[3]*r0s2mr2s0, +q[0]*r1s3mr3s1-q[1]*r0s3mr3s0+q[3]*r0s1mr1s0, -q[0]*r1s2mr2s1+q[1]*r0s2mr2s0-q[2]*r0s1mr1s0); } template inline double CosPhi(const Vec4 &pi,const Vec4 &pj, const Vec4 &pk,const Vec4 &pl) { Scalar sij((pi+pj).Abs2()), sik((pi+pk).Abs2()), sil((pi+pl).Abs2()); Scalar sjk((pj+pk).Abs2()), sjl((pj+pl).Abs2()), skl((pk+pl).Abs2()); Scalar si(pi.Abs2()), sj(pj.Abs2()), sk(pk.Abs2()), sl(pl.Abs2()); Scalar cp(skl*(sik*sjl+sil*sjk)-sk*sil*sjl-sl*sik*sjk-sij*skl*skl+sij*sk*sl); return cp/sqrt((2.0*skl*sik*sil-sk*sil*sil-sl*sik*sik-si*skl*skl+si*sk*sl)* (2.0*skl*sjk*sjl-sk*sjl*sjl-sl*sjk*sjk-sj*skl*skl+sj*sk*sl)); } template std::ostream& operator<<(std::ostream& s, const Vec4& vec) { return s<<'('< \brief implementation of a 4 dimensional Minkowski vector and its algebraic operations This class can be used as Minkowski vector with arbitrary scalar types. All necessary operations, e.g. addition, scalar product, cross product, etc. are available. If you want to use this vector for a new scalar type, you might have to implement specific functions for this type in MathTools.H. "Fallback" types for operations with two vectors of different types, e. g. "complex and double become complex" have to be specified in MathTools.H as well, see the DECLARE_PROMOTE macro. */ /*! \fn Vec4::Vec4() \brief Standard Constructor */ /*! \fn Vec4::Vec4(Scalar x0,Scalar x1, Scalar x2, Scalar x3){ \brief Special Constructor, templated in Scalar, taking 4 single components. */ /*! \fn Vec4::Vec4(const Scalar E, const Vec3D& v) \brief Special Constructor taking the space part and a Energie */ /*! \fn inline const Scalar Vec4::Abs2() const \brief returns \f$ x_0^2 - (x_1^2 + x_2^2 + x_3^2) \f$. */ /*! \fn inline Scalar& Vec4::operator[] (int i) \brief returns \f$x_i\f$. (May be manipulated.) */ /*! \fn inline const Scalar Vec4::operator[] (int i) const \brief returns \f$x_i\f$. */ #endif