Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2025-12-15 10:30:40 

0001 #ifndef ATOOLS_Math_Vec4_H
0002 #define ATOOLS_Math_Vec4_H
0003 
0004 #include "ATOOLS/Math/MathTools.H"
0005 
0006 namespace ATOOLS {
0007   template<typename Scalar> class Vec4;
0008   template<typename Scalar> class Vec3;
0009 
0010   template<typename Scalar>
0011   class Vec4 {
0012     Scalar m_x[4];
0013     template <typename Scalar2> friend class Vec4;
0014   public:
0015     inline Vec4() {
0016       m_x[0]=m_x[1]=m_x[2]=m_x[3]=Scalar(0.0);
0017     }
0018     template<typename Scalar2>
0019     inline Vec4(const Vec4<Scalar2> & vec) {
0020       m_x[0] = vec[0]; m_x[1] = vec[1];
0021       m_x[2] = vec[2]; m_x[3] = vec[3];
0022     }
0023     /*!
0024      \brief Special Constructor, templated in Scalar, taking 4 single components.
0025     */
0026     inline Vec4(const Scalar& x0, const Scalar& x1,
0027                 const Scalar& x2, const Scalar& x3) {
0028       m_x[0] = x0; m_x[1] = x1; m_x[2] = x2; m_x[3] = x3;
0029     }
0030     /*!
0031      \brief Special Constructor taking the space part and a Energie
0032     */
0033     inline Vec4(const Scalar& E, const Vec3<Scalar>& xn) {
0034       m_x[0] = E;
0035       m_x[1] = xn[1];
0036       m_x[2] = xn[2];
0037       m_x[3] = xn[3];
0038     }
0039 
0040     inline Scalar& operator[] (int i) {
0041       return m_x[i];
0042     }
0043 
0044     inline const Scalar operator[] (int i) const {
0045       return m_x[i];
0046     }
0047 
0048     inline Vec4<Scalar>& operator+=(const Vec4<Scalar>& v) {
0049       m_x[0] += v[0];
0050       m_x[1] += v[1];
0051       m_x[2] += v[2];
0052       m_x[3] += v[3];
0053       return *this;
0054     }
0055 
0056     inline Vec4<Scalar>& operator-=(const Vec4<Scalar>& v) {
0057       m_x[0] -= v[0];
0058       m_x[1] -= v[1];
0059       m_x[2] -= v[2];
0060       m_x[3] -= v[3];
0061       return *this;
0062     }
0063 
0064     inline Vec4<Scalar>& operator*= (const Scalar& scal) {
0065       m_x[0] *= scal;
0066       m_x[1] *= scal;
0067       m_x[2] *= scal;
0068       m_x[3] *= scal;
0069       return *this;
0070     }
0071 
0072     inline Vec4<Scalar> operator-() const {
0073       return Vec4<Scalar>(-m_x[0],-m_x[1],-m_x[2],-m_x[3]);
0074     }
0075 
0076     template<typename Scalar2> inline PROMOTE(Scalar,Scalar2)
0077     operator*(const Vec4<Scalar2>& v2) const {
0078       return m_x[0]*v2[0]-m_x[1]*v2[1]-m_x[2]*v2[2]-m_x[3]*v2[3];
0079     }
0080 
0081     template<typename Scalar2> inline Vec4<PROMOTE(Scalar,Scalar2)>
0082     operator+ (const Vec4<Scalar2>& v2) const {
0083       return Vec4<PROMOTE(Scalar,Scalar2)>
0084         (m_x[0]+v2[0],m_x[1]+v2[1],m_x[2]+v2[2],m_x[3]+v2[3]);
0085     }
0086 
0087     template<typename Scalar2> inline Vec4<PROMOTE(Scalar,Scalar2)>
0088     operator- (const Vec4<Scalar2>& v2) const {
0089       return Vec4<PROMOTE(Scalar,Scalar2)>
0090         (m_x[0]-v2[0],m_x[1]-v2[1],m_x[2]-v2[2],m_x[3]-v2[3]);
0091     }
0092 
0093     template<typename Scalar2> inline Vec4<PROMOTE(Scalar,Scalar2)>
0094     operator* (const Scalar2& s) const {
0095       return Vec4<PROMOTE(Scalar,Scalar2)>
0096         (s*m_x[0],s*m_x[1],s*m_x[2],s*m_x[3]);
0097     }
0098 
0099     template<typename Scalar2> inline Vec4<PROMOTE(Scalar,Scalar2)>
0100     operator/ (const Scalar2& scal) const {
0101       return (*this)*(1./scal);
0102     }
0103 
0104 
0105     inline bool Nan() const {
0106       for(short unsigned int i(0);i<4;++i)
0107         if (ATOOLS::IsNan<Scalar>(m_x[i])) return true;
0108       return false;
0109     }
0110     inline bool IsZero() const {
0111       for(short unsigned int i(0);i<4;++i) 
0112         if (!ATOOLS::IsZero<Scalar>(m_x[i])) return false;
0113       return true;
0114     }
0115 
0116 
0117     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]; }      
0118     Scalar RelAbs2() const {
0119       const auto abs2 = Abs2();
0120       return (abs2 == 0) ? 0 : abs2/(m_x[0]*m_x[0]);
0121     }
0122     Scalar Abs() const { return sqrt(Abs2()); }
0123     double Mass() const {
0124       return sqrt(ATOOLS::Abs<Scalar>(Abs2()));
0125     }
0126     Scalar P() const { return PSpat(); }
0127 
0128     Scalar PPerp2() const { return m_x[1]*m_x[1]+m_x[2]*m_x[2]; }
0129     Scalar PPerp() const { return sqrt(PPerp2()); }
0130     Scalar MPerp2() const { return (m_x[0]+m_x[3])*(m_x[0]-m_x[3]); }
0131     Scalar MPerp() const { return sqrt(MPerp2()); }
0132     Scalar MPerp2(const Vec4<Scalar> &ref) const {
0133       return Abs2()+PPerp2(ref);
0134     }
0135     Scalar MPerp(const Vec4<Scalar> &ref) const {
0136       return sqrt(MPerp2(ref));
0137     }
0138     Scalar EPerp2() const {
0139       return m_x[0]*m_x[0]*PPerp2()/PSpat2();
0140     }
0141     Scalar EPerp() const { return sqrt(EPerp2()); }
0142     Scalar PPlus() const { return m_x[0]+m_x[3]; }
0143     Scalar PMinus() const { return m_x[0]-m_x[3]; }
0144     Scalar E() const { return m_x[0]; }
0145 
0146     Scalar PSpat2() const {
0147       return m_x[1]*m_x[1]+m_x[2]*m_x[2]+m_x[3]*m_x[3];
0148     }
0149     Scalar PSpat() const { return sqrt(PSpat2()); }
0150     Scalar Y() const { return 0.5*log(PPlus()/PMinus()); }
0151 
0152     inline Vec4<Scalar> Perp() const { 
0153       return Vec4<Scalar>(0.,m_x[1],m_x[2],0.);
0154     }
0155     inline Vec4<Scalar> Long() const { 
0156       return Vec4<Scalar>(m_x[0],0.,0.,m_x[3]);
0157     }
0158     inline Vec4<Scalar> Plus() const {
0159       Scalar pplus=0.5*PPlus();
0160       return Vec4<Scalar>(pplus,0.0,0.0,pplus);
0161     }
0162     inline Vec4<Scalar> Minus() const {
0163       Scalar pminus=0.5*PMinus();
0164       return Vec4<Scalar>(pminus,0.0,0.0,-pminus);
0165     }
0166     
0167     Scalar PPerp2(const Vec4<Scalar>& ref) const {
0168       Scalar ppref = PPerp(ref);
0169       return ppref*ppref;
0170     }
0171     Scalar PPerp(const Vec4<Scalar>& ref) const {
0172       Vec3<Scalar> perp=1./ATOOLS::Max(Vec3<Scalar>(ref).Abs(),1.e-12)*
0173         Vec3<Scalar>(ref);
0174       perp=Vec3<Scalar>(*this)-perp*(perp*Vec3<Scalar>(*this));
0175       return perp.Abs();
0176     }
0177 
0178     // some functions which only make sense for doubles
0179     // and thus have to be specialised for the type
0180     Scalar CosPhi() const;
0181     Scalar SinPhi() const;
0182     Scalar Phi() const;
0183     Scalar CosTheta() const;
0184     Scalar SinTheta() const;
0185     Scalar Theta() const;
0186     Scalar Eta() const;
0187     Scalar CosTheta(const Vec4<Scalar>& ref) const;
0188     Scalar Theta(const Vec4<Scalar>& ref) const;
0189     Scalar Eta(const Vec4<Scalar>& ref) const;
0190     Scalar CosDPhi(const Vec4<Scalar>& ref) const;
0191     Scalar DPhi(const Vec4<Scalar>& ref) const;
0192     Scalar DEta(const Vec4<Scalar>& ref) const;
0193     Scalar DY(const Vec4<Scalar>& ref) const;
0194     Scalar DR(const Vec4<Scalar>& ref) const;
0195     Scalar DR2(const Vec4<Scalar>& ref) const;
0196     Scalar DRy(const Vec4<Scalar>& ref) const;
0197     Scalar DR2y(const Vec4<Scalar>& ref) const;
0198     Scalar SmallOMCT(const Vec4<Scalar>& ref) const;
0199     Scalar SmallMLDP(const Vec4<Scalar>& ref) const;
0200 
0201     // standard vectors:
0202     const static Vec4<Scalar> XVEC;
0203     const static Vec4<Scalar> YVEC;
0204     const static Vec4<Scalar> ZVEC;
0205   };
0206 
0207   template<typename Scalar,typename Scal2> inline Vec4<PROMOTE(Scalar,Scal2)>
0208   operator* (const Scal2& s, const Vec4<Scalar>& v) {
0209     return Vec4<PROMOTE(Scalar,Scal2)>(v*s);
0210   }
0211 
0212   template<typename S1,typename S2,typename S3> inline
0213   Vec4<PROMOTE(S1,PROMOTE(S2,S3))>
0214   cross(const Vec4<S1>&q, const Vec4<S2>&r, const Vec4<S3>&s)
0215   {
0216     PROMOTE(S1,PROMOTE(S2,S3)) r0s1mr1s0(r[0]*s[1]-r[1]*s[0]), 
0217       r0s2mr2s0(r[0]*s[2]-r[2]*s[0]), r0s3mr3s0(r[0]*s[3]-r[3]*s[0]),
0218       r1s2mr2s1(r[1]*s[2]-r[2]*s[1]), r1s3mr3s1(r[1]*s[3]-r[3]*s[1]),
0219       r2s3mr3s2(r[2]*s[3]-r[3]*s[2]);
0220     return Vec4<PROMOTE(S1,PROMOTE(S2,S3))>
0221       (-q[1]*r2s3mr3s2+q[2]*r1s3mr3s1-q[3]*r1s2mr2s1,
0222        -q[0]*r2s3mr3s2+q[2]*r0s3mr3s0-q[3]*r0s2mr2s0,
0223        +q[0]*r1s3mr3s1-q[1]*r0s3mr3s0+q[3]*r0s1mr1s0,
0224        -q[0]*r1s2mr2s1+q[1]*r0s2mr2s0-q[2]*r0s1mr1s0);
0225   }
0226 
0227   template<typename Scalar> inline
0228   double CosPhi(const Vec4<Scalar> &pi,const Vec4<Scalar> &pj,
0229         const Vec4<Scalar> &pk,const Vec4<Scalar> &pl)
0230   {
0231     Scalar sij((pi+pj).Abs2()), sik((pi+pk).Abs2()), sil((pi+pl).Abs2());
0232     Scalar sjk((pj+pk).Abs2()), sjl((pj+pl).Abs2()), skl((pk+pl).Abs2());
0233     Scalar si(pi.Abs2()), sj(pj.Abs2()), sk(pk.Abs2()), sl(pl.Abs2());
0234     Scalar cp(skl*(sik*sjl+sil*sjk)-sk*sil*sjl-sl*sik*sjk-sij*skl*skl+sij*sk*sl);
0235     return cp/sqrt((2.0*skl*sik*sil-sk*sil*sil-sl*sik*sik-si*skl*skl+si*sk*sl)*
0236            (2.0*skl*sjk*sjl-sk*sjl*sjl-sl*sjk*sjk-sj*skl*skl+sj*sk*sl));
0237   }
0238   
0239   template<typename Scalar>
0240   std::ostream& operator<<(std::ostream& s, const Vec4<Scalar>& vec) {
0241     return s<<'('<<vec[0]<<','<<vec[1]<<','<<vec[2]<<','<<vec[3]<<')';
0242   }
0243 
0244 }
0245 
0246 /*!
0247  \file
0248  \brief   contains class Vec4
0249 */
0250 
0251 /*!
0252  \class Vec4<Scalar>
0253  \brief implementation of a 4 dimensional Minkowski vector and its algebraic
0254  operations
0255 
0256  This class can be used as Minkowski vector with arbitrary scalar types.
0257  All necessary operations, e.g. addition, scalar product, cross product, 
0258  etc. are available.
0259  If you want to use this vector for a new scalar type, you might have to
0260  implement specific functions for this type in MathTools.H.
0261  "Fallback" types for operations with two vectors of different types, e. g.
0262  "complex and double become complex" have to be specified in MathTools.H as
0263  well, see the DECLARE_PROMOTE macro.
0264 */
0265 
0266 /*!
0267  \fn Vec4::Vec4()
0268  \brief Standard Constructor
0269 */
0270 
0271 /*!
0272  \fn Vec4::Vec4(const Scalar &x0, const Scalar &x1, const Scalar &x2, const Scalar &x3){
0273  \brief Special Constructor, templated in Scalar, taking 4 single components.
0274 */
0275 
0276 /*!
0277  \fn ATOOLS::Vec4<Scalar>::Vec4(const Scalar &E, const Vec3<Scalar>& v)
0278  \brief Special Constructor taking the space part and a Energie
0279 */
0280 
0281 /*!
0282  \fn inline const Scalar Vec4::Abs2() const
0283  \brief returns \f$ x_0^2 - (x_1^2 + x_2^2 + x_3^2) \f$.
0284 */
0285 
0286 /*!
0287  \fn   inline Scalar& Vec4::operator[] (int i)
0288  \brief returns \f$x_i\f$. (May be manipulated.)
0289 */
0290 
0291 /*!
0292  \fn   inline const Scalar Vec4::operator[] (int i) const
0293  \brief returns \f$x_i\f$.
0294 */
0295 
0296 #endif
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!