Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2026-07-26 09:24:23

0001 // @(#)root/mathcore:$Id$
0002 // Authors: W. Brown, M. Fischler, L. Moneta    2005
0003 
0004  /**********************************************************************
0005   *                                                                    *
0006   * Copyright (c) 2005 , LCG ROOT MathLib Team                         *
0007   *                                                                    *
0008   *                                                                    *
0009   **********************************************************************/
0010 
0011 // Header file for class PositionVector3D
0012 //
0013 // Created by: Lorenzo Moneta  at Mon May 30 15:25:04 2005
0014 //
0015 // Last update: $Id$
0016 //
0017 #ifndef ROOT_Math_GenVector_PositionVector3D
0018 #define ROOT_Math_GenVector_PositionVector3D  1
0019 
0020 #include "Math/GenVector/DisplacementVector3Dfwd.h"
0021 
0022 #include "Math/GenVector/Cartesian3D.h"
0023 
0024 #include "Math/GenVector/GenVectorIO.h"
0025 
0026 #include "Math/GenVector/BitReproducible.h"
0027 
0028 #include "Math/GenVector/CoordinateSystemTags.h"
0029 
0030 
0031 #include <cassert>
0032 
0033 namespace ROOT {
0034 
0035   namespace Math {
0036 
0037 
0038 //__________________________________________________________________________________________
0039     /**
0040      Class describing a generic position vector (point) in 3 dimensions.
0041      This class is templated on the type of Coordinate system.
0042      One example is the XYZPoint which is a vector based on
0043      double precision x,y,z data members by using the
0044      ROOT::Math::Cartesian3D<double> Coordinate system.
0045      The class is having also an extra template parameter, the coordinate system tag,
0046      to be able to identify (tag) vector described in different reference coordinate system,
0047      like global or local coordinate systems.
0048 
0049      @ingroup GenVector
0050 
0051      @see GenVector
0052     */
0053 
0054     template <class CoordSystem, class Tag = DefaultCoordinateSystemTag >
0055     class PositionVector3D {
0056 
0057     public:
0058 
0059       typedef typename CoordSystem::Scalar Scalar;
0060       typedef CoordSystem CoordinateType;
0061       typedef Tag  CoordinateSystemTag;
0062 
0063       // ------ ctors ------
0064 
0065       /**
0066          Default constructor. Construct an empty object with zero values
0067       */
0068 
0069       constexpr PositionVector3D() noexcept = default;
0070 
0071       /**
0072          Construct from three values of type <em>Scalar</em>.
0073          In the case of a XYZPoint the values are x,y,z
0074          In the case of  a polar vector they are r,theta,phi
0075       */
0076       constexpr PositionVector3D(const Scalar &a, const Scalar &b, const Scalar &c) noexcept : fCoordinates(a, b, c) {}
0077 
0078      /**
0079           Construct from a position vector expressed in different
0080           coordinates, or using a different Scalar type
0081       */
0082       template <class T>
0083       explicit constexpr PositionVector3D( const PositionVector3D<T,Tag> & v) :
0084         fCoordinates ( v.Coordinates() ) { }
0085 
0086      /**
0087           Construct from an arbitrary displacement vector
0088       */
0089       template <class T>
0090       explicit constexpr PositionVector3D( const DisplacementVector3D<T,Tag> & p) :
0091         fCoordinates ( p.Coordinates() ) { }
0092 
0093       /**
0094           Construct from a foreign 3D vector type, for example, Hep3Vector
0095           Precondition: v must implement methods x(), y() and z()
0096       */
0097       template <class ForeignVector>
0098       explicit constexpr PositionVector3D( const ForeignVector & v) :
0099         fCoordinates ( Cartesian3D<Scalar>( v.x(), v.y(), v.z() ) ) { }
0100 
0101 #ifdef LATER
0102       /**
0103          construct from a generic linear algebra  vector of at least size 3
0104          implementing operator []. This could be also a C array
0105          \par v  LAVector
0106          \par index0   index where coordinates starts (typically zero)
0107          It works for all Coordinates types,
0108          ( x= v[index0] for Cartesian and r=v[index0] for Polar )
0109       */
0110       template <class LAVector>
0111       PositionVector3D(const LAVector & v, size_t index0 ) {
0112         fCoordinates = CoordSystem  ( v[index0], v[index0+1], v[index0+2] );
0113       }
0114 #endif
0115 
0116       // compiler-generated copy ctor and dtor are fine.
0117 
0118       // ------ assignment ------
0119 
0120       /**
0121           Assignment operator from a position vector of arbitrary type
0122       */
0123       template <class OtherCoords>
0124       PositionVector3D & operator=
0125                         ( const PositionVector3D<OtherCoords,Tag> & v) {
0126         fCoordinates = v.Coordinates();
0127         return *this;
0128       }
0129 
0130       /**
0131           Assignment operator from a displacement vector of arbitrary type
0132       */
0133       template <class OtherCoords>
0134       PositionVector3D & operator=
0135                         ( const DisplacementVector3D<OtherCoords,Tag> & v) {
0136         fCoordinates = v.Coordinates();
0137         return *this;
0138       }
0139 
0140       /**
0141           Assignment from a foreign 3D vector type, for example, Hep3Vector
0142           Precondition: v must implement methods x(), y() and z()
0143       */
0144       template <class ForeignVector>
0145       PositionVector3D & operator= ( const ForeignVector & v) {
0146         SetXYZ( v.x(),  v.y(), v.z() );
0147         return *this;
0148       }
0149 
0150 #ifdef LATER
0151       /**
0152          assign from a generic linear algebra  vector of at least size 3
0153          implementing operator [].
0154          \par v  LAVector
0155          \par index0   index where coordinates starts (typically zero)
0156          It works for all Coordinates types,
0157          ( x= v[index0] for Cartesian and r=v[index0] for Polar )
0158       */
0159       template <class LAVector>
0160       PositionVector3D & assignFrom(const LAVector & v, size_t index0 = 0) {
0161         fCoordinates = CoordSystem  ( v[index0], v[index0+1], v[index0+2] );
0162         return *this;
0163       }
0164 #endif
0165 
0166       /**
0167           Retrieve a copy of the coordinates object
0168       */
0169       const CoordSystem & Coordinates() const {
0170         return fCoordinates;
0171       }
0172 
0173       /**
0174          Set internal data based on a C-style array of 3 Scalar numbers
0175        */
0176       PositionVector3D<CoordSystem, Tag>& SetCoordinates( const Scalar src[] )
0177        { fCoordinates.SetCoordinates(src); return *this;  }
0178 
0179       /**
0180          Set internal data based on 3 Scalar numbers
0181        */
0182       PositionVector3D<CoordSystem, Tag>& SetCoordinates( Scalar a, Scalar b, Scalar c )
0183        { fCoordinates.SetCoordinates(a, b, c); return *this; }
0184 
0185       /**
0186          Set internal data based on 3 Scalars at *begin to *end
0187        */
0188       template <class IT>
0189       PositionVector3D<CoordSystem, Tag>& SetCoordinates( IT begin, IT end )
0190       { IT a = begin; IT b = ++begin; IT c = ++begin;
0191         (void)end;
0192         assert (++begin==end);
0193         SetCoordinates (*a,*b,*c);
0194         return *this;
0195       }
0196 
0197       /**
0198         get internal data into 3 Scalar numbers
0199        */
0200       void GetCoordinates( Scalar& a, Scalar& b, Scalar& c ) const
0201                             { fCoordinates.GetCoordinates(a, b, c);  }
0202 
0203       /**
0204          get internal data into a C-style array of 3 Scalar numbers
0205        */
0206       void GetCoordinates( Scalar dest[] ) const
0207                             { fCoordinates.GetCoordinates(dest);  }
0208 
0209       /**
0210          get internal data into 3 Scalars at *begin to *end (3 past begin)
0211        */
0212       template <class IT>
0213       void GetCoordinates( IT begin, IT end ) const
0214       { IT a = begin; IT b = ++begin; IT c = ++begin;
0215         (void)end;
0216         assert (++begin==end);
0217         GetCoordinates (*a,*b,*c);
0218       }
0219 
0220       /**
0221          get internal data into 3 Scalars at *begin
0222        */
0223       template <class IT>
0224       void GetCoordinates( IT begin ) const {
0225          Scalar a = Scalar(0);
0226          Scalar b = Scalar(0);
0227          Scalar c = Scalar(0);
0228          GetCoordinates(a, b, c);
0229          *begin++ = a;
0230          *begin++ = b;
0231          *begin   = c;
0232       }
0233 
0234       /**
0235          set the values of the vector from the cartesian components (x,y,z)
0236          (if the vector is held in polar or cylindrical eta coordinates,
0237          then (x, y, z) are converted to that form)
0238        */
0239       PositionVector3D<CoordSystem, Tag>& SetXYZ (Scalar a, Scalar b, Scalar c) {
0240             fCoordinates.SetXYZ(a,b,c);
0241             return *this;
0242       }
0243 
0244       // ------------------- Equality -----------------
0245 
0246       /**
0247         Exact equality
0248        */
0249       bool operator==(const PositionVector3D & rhs) const {
0250         return fCoordinates==rhs.fCoordinates;
0251       }
0252       bool operator!= (const PositionVector3D & rhs) const {
0253         return !(operator==(rhs));
0254       }
0255 
0256       // ------ Individual element access, in various coordinate systems ------
0257 
0258       /**
0259           Dimension
0260       */
0261       unsigned int Dimension() const
0262       {
0263          return fDimension;
0264       };
0265 
0266       /**
0267           Cartesian X, converting if necessary from internal coordinate system.
0268       */
0269       Scalar X() const { return fCoordinates.X(); }
0270 
0271       /**
0272           Cartesian Y, converting if necessary from internal coordinate system.
0273       */
0274       Scalar Y() const { return fCoordinates.Y(); }
0275 
0276       /**
0277           Cartesian Z, converting if necessary from internal coordinate system.
0278       */
0279       Scalar Z() const { return fCoordinates.Z(); }
0280 
0281       /**
0282           Polar R, converting if necessary from internal coordinate system.
0283       */
0284       Scalar R() const { return fCoordinates.R(); }
0285 
0286       /**
0287           Polar theta, converting if necessary from internal coordinate system.
0288       */
0289       Scalar Theta() const { return fCoordinates.Theta(); }
0290 
0291       /**
0292           Polar phi, converting if necessary from internal coordinate system.
0293       */
0294       Scalar Phi() const { return fCoordinates.Phi(); }
0295 
0296       /**
0297           Polar eta, converting if necessary from internal coordinate system.
0298       */
0299       Scalar Eta() const { return fCoordinates.Eta(); }
0300 
0301       /**
0302           Cylindrical transverse component rho
0303       */
0304       Scalar Rho() const { return fCoordinates.Rho(); }
0305 
0306       // ----- Other fundamental properties -----
0307 
0308       /**
0309           Magnitute squared ( r^2 in spherical coordinate)
0310       */
0311       Scalar Mag2() const { return fCoordinates.Mag2();}
0312 
0313       /**
0314          Transverse component squared (rho^2 in cylindrical coordinates.
0315       */
0316       Scalar Perp2() const { return fCoordinates.Perp2();}
0317 
0318       // It is physically meaningless to speak of the unit vector corresponding
0319       // to a point.
0320 
0321       // ------ Setting individual elements present in coordinate system ------
0322 
0323       /**
0324          Change X - Cartesian3D coordinates only
0325       */
0326        PositionVector3D<CoordSystem, Tag>& SetX (Scalar xx) { fCoordinates.SetX(xx); return *this;}
0327 
0328       /**
0329          Change Y - Cartesian3D coordinates only
0330       */
0331        PositionVector3D<CoordSystem, Tag>& SetY (Scalar yy) { fCoordinates.SetY(yy); return *this;}
0332 
0333       /**
0334          Change Z - Cartesian3D coordinates only
0335       */
0336        PositionVector3D<CoordSystem, Tag>& SetZ (Scalar zz) { fCoordinates.SetZ(zz); return *this;}
0337 
0338       /**
0339          Change R - Polar3D coordinates only
0340       */
0341        PositionVector3D<CoordSystem, Tag>& SetR (Scalar rr) { fCoordinates.SetR(rr); return *this;}
0342 
0343       /**
0344          Change Theta - Polar3D coordinates only
0345       */
0346        PositionVector3D<CoordSystem, Tag>& SetTheta (Scalar ang) { fCoordinates.SetTheta(ang); return *this;}
0347 
0348       /**
0349          Change Phi - Polar3D or CylindricalEta3D coordinates
0350       */
0351        PositionVector3D<CoordSystem, Tag>& SetPhi (Scalar ang) { fCoordinates.SetPhi(ang); return *this;}
0352 
0353       /**
0354          Change Rho - CylindricalEta3D coordinates only
0355       */
0356        PositionVector3D<CoordSystem, Tag>& SetRho (Scalar rr) { fCoordinates.SetRho(rr); return *this;}
0357 
0358       /**
0359          Change Eta - CylindricalEta3D coordinates only
0360       */
0361        PositionVector3D<CoordSystem, Tag>& SetEta (Scalar etaval) { fCoordinates.SetEta(etaval); return *this;}
0362 
0363       // ------ Operations combining two vectors ------
0364       // need to specialize to exclude those with a different tags
0365 
0366      /**
0367       Return the scalar (Dot) product of this with a displacement vector in
0368       any coordinate system, but with the same tag
0369       */
0370       template< class OtherCoords >
0371       Scalar Dot( const  DisplacementVector3D<OtherCoords,Tag> & v) const {
0372         return X()*v.x() + Y()*v.y() + Z()*v.z();
0373       }
0374 
0375 
0376       /**
0377          Return vector (Cross) product of this point with a displacement, as a
0378          point vector in this coordinate system of the first.
0379       */
0380       template< class OtherCoords >
0381       PositionVector3D Cross( const DisplacementVector3D<OtherCoords,Tag> & v) const  {
0382         PositionVector3D  result;
0383         result.SetXYZ (  Y()*v.z() - v.y()*Z(),
0384                          Z()*v.x() - v.z()*X(),
0385                          X()*v.y() - v.x()*Y() );
0386         return result;
0387       }
0388 
0389       // The Dot and Cross products of a pair of point vectors are physically
0390       // meaningless concepts and thus are defined as private methods
0391 
0392       // It is physically meaningless to speak of the Unit vector corresponding
0393       // to a point.
0394 
0395 
0396       /**
0397           Self Addition with a displacement vector.
0398       */
0399       template <class OtherCoords>
0400       PositionVector3D & operator+= (const  DisplacementVector3D<OtherCoords,Tag> & v)
0401       {
0402         SetXYZ( X() + v.X(), Y() + v.Y(), Z() + v.Z() );
0403         return *this;
0404       }
0405 
0406       /**
0407           Self Difference with a displacement vector.
0408       */
0409       template <class OtherCoords>
0410       PositionVector3D & operator-= (const  DisplacementVector3D<OtherCoords,Tag> & v)
0411       {
0412         SetXYZ(  X() - v.X(), Y() - v.Y(), Z() - v.Z() );
0413         return *this;
0414       }
0415 
0416       /**
0417          multiply this vector by a scalar quantity
0418       */
0419       PositionVector3D & operator *= (Scalar a) {
0420         fCoordinates.Scale(a);
0421         return *this;
0422       }
0423 
0424       /**
0425          divide this vector by a scalar quantity
0426       */
0427       PositionVector3D & operator /= (Scalar a) {
0428         fCoordinates.Scale(1/a);
0429         return *this;
0430       }
0431 
0432       // The following methods (v*a and v/a) could instead be free functions.
0433       // They were moved into the class to solve a problem on AIX.
0434       /**
0435         Multiply a vector by a real number
0436       */
0437       PositionVector3D operator * ( Scalar a ) const {
0438         PositionVector3D tmp(*this);
0439         tmp *= a;
0440         return tmp;
0441       }
0442 
0443       /**
0444          Division of a vector with a real number
0445        */
0446       PositionVector3D operator / (Scalar a) const {
0447         PositionVector3D tmp(*this);
0448         tmp /= a;
0449         return tmp;
0450       }
0451 
0452       // Limited backward name compatibility with CLHEP
0453 
0454       Scalar x()     const { return fCoordinates.X();     }
0455       Scalar y()     const { return fCoordinates.Y();     }
0456       Scalar z()     const { return fCoordinates.Z();     }
0457       Scalar r()     const { return fCoordinates.R();     }
0458       Scalar theta() const { return fCoordinates.Theta(); }
0459       Scalar phi()   const { return fCoordinates.Phi();   }
0460       Scalar eta()   const { return fCoordinates.Eta();   }
0461       Scalar rho()   const { return fCoordinates.Rho();   }
0462       Scalar mag2()  const { return fCoordinates.Mag2();  }
0463       Scalar perp2() const { return fCoordinates.Perp2(); }
0464 
0465     private:
0466 
0467       CoordSystem fCoordinates;
0468       static constexpr unsigned int fDimension = CoordinateType::Dimension;
0469 
0470       // Prohibited methods
0471 
0472       // this should not compile (if from a vector or points with different tag
0473 
0474       template <class OtherCoords, class OtherTag>
0475       explicit constexpr PositionVector3D( const PositionVector3D<OtherCoords, OtherTag> & );
0476 
0477       template <class OtherCoords, class OtherTag>
0478       explicit constexpr PositionVector3D( const DisplacementVector3D<OtherCoords, OtherTag> & );
0479 
0480       template <class OtherCoords, class OtherTag>
0481       PositionVector3D & operator=( const PositionVector3D<OtherCoords, OtherTag> & );
0482 
0483       template <class OtherCoords, class OtherTag>
0484       PositionVector3D & operator=( const DisplacementVector3D<OtherCoords, OtherTag> & );
0485 
0486       template <class OtherCoords, class OtherTag>
0487       PositionVector3D & operator+=(const  DisplacementVector3D<OtherCoords, OtherTag> & );
0488 
0489       template <class OtherCoords, class OtherTag>
0490       PositionVector3D & operator-=(const  DisplacementVector3D<OtherCoords, OtherTag> & );
0491 
0492 //       /**
0493 //          Dot product of two position vectors is inappropriate
0494 //       */
0495 //       template <class T2, class U>
0496 //       PositionVector3D Dot( const PositionVector3D<T2,U> & v) const;
0497 
0498 //       /**
0499 //          Cross product of two position vectors is inappropriate
0500 //       */
0501 //       template <class T2, class U>
0502 //       PositionVector3D Cross( const PositionVector3D<T2,U> & v) const;
0503 
0504 
0505 
0506     };
0507 
0508 // ---------- PositionVector3D class template ends here ----------------
0509 // ---------------------------------------------------------------------
0510 
0511     /**
0512        Multiplication of a position vector by real number  a*v
0513     */
0514     template <class CoordSystem, class U>
0515     inline
0516     PositionVector3D<CoordSystem>
0517     operator * ( typename PositionVector3D<CoordSystem,U>::Scalar a,
0518                  PositionVector3D<CoordSystem,U> v) {
0519       return v *= a;
0520       // Note - passing v by value and using operator *= may save one
0521       // copy relative to passing v by const ref and creating a temporary.
0522     }
0523 
0524     /**
0525         Difference between two PositionVector3D vectors.
0526         The result is a DisplacementVector3D.
0527         The (coordinate system) type of the returned vector is defined to
0528         be identical to that of the first position vector.
0529     */
0530 
0531     template <class CoordSystem1, class CoordSystem2, class U>
0532     inline
0533     DisplacementVector3D<CoordSystem1,U>
0534     operator-( const PositionVector3D<CoordSystem1,U> & v1,
0535                const PositionVector3D<CoordSystem2,U> & v2) {
0536       return DisplacementVector3D<CoordSystem1,U>( Cartesian3D<typename CoordSystem1::Scalar>(
0537                                                                                v1.X()-v2.X(), v1.Y()-v2.Y(),v1.Z()-v2.Z() )
0538                                              );
0539     }
0540 
0541     /**
0542         Addition of a PositionVector3D and a DisplacementVector3D.
0543         The return type is a PositionVector3D,
0544         of the same (coordinate system) type as the input PositionVector3D.
0545     */
0546     template <class CoordSystem1, class CoordSystem2, class U>
0547     inline
0548     PositionVector3D<CoordSystem2,U>
0549     operator+( PositionVector3D<CoordSystem2,U> p1,
0550                const DisplacementVector3D<CoordSystem1,U>  & v2)        {
0551       return p1 += v2;
0552     }
0553 
0554     /**
0555         Addition of a DisplacementVector3D and a PositionVector3D.
0556         The return type is a PositionVector3D,
0557         of the same (coordinate system) type as the input PositionVector3D.
0558     */
0559     template <class CoordSystem1, class CoordSystem2, class U>
0560     inline
0561     PositionVector3D<CoordSystem2,U>
0562     operator+( DisplacementVector3D<CoordSystem1,U> const & v1,
0563                PositionVector3D<CoordSystem2,U> p2)        {
0564       return p2 += v1;
0565     }
0566 
0567     /**
0568         Subtraction of a DisplacementVector3D from a PositionVector3D.
0569         The return type is a PositionVector3D,
0570         of the same (coordinate system) type as the input PositionVector3D.
0571     */
0572     template <class CoordSystem1, class CoordSystem2, class U>
0573     inline
0574     PositionVector3D<CoordSystem2,U>
0575     operator-( PositionVector3D<CoordSystem2,U> p1,
0576                DisplacementVector3D<CoordSystem1,U> const & v2)        {
0577       return p1 -= v2;
0578     }
0579 
0580     // Scaling of a position vector with a real number is not physically meaningful
0581 
0582     // ------------- I/O to/from streams -------------
0583 
0584     template <
0585        class char_t, class traits_t, class T, class U,
0586        typename std::enable_if<std::is_arithmetic<typename PositionVector3D<T, U>::Scalar>::value>::type * = nullptr>
0587     std::basic_ostream<char_t, traits_t> &operator<<(std::basic_ostream<char_t, traits_t> &os,
0588                                                      PositionVector3D<T, U> const &v)
0589     {
0590        if (os) {
0591 
0592           typename T::Scalar a = 0;
0593           typename T::Scalar b = 0;
0594           typename T::Scalar c = 0;
0595           v.GetCoordinates(a, b, c);
0596 
0597           if (detail::get_manip(os, detail::bitforbit)) {
0598              detail::set_manip(os, detail::bitforbit, '\00');
0599              typedef GenVector_detail::BitReproducible BR;
0600              BR::Output(os, a);
0601              BR::Output(os, b);
0602              BR::Output(os, c);
0603           } else {
0604              os << detail::get_manip(os, detail::open) << a << detail::get_manip(os, detail::sep) << b
0605                 << detail::get_manip(os, detail::sep) << c << detail::get_manip(os, detail::close);
0606           }
0607       }
0608       return os;
0609     }  // op<< <>()
0610 
0611     template <
0612        class char_t, class traits_t, class T, class U,
0613        typename std::enable_if<!std::is_arithmetic<typename PositionVector3D<T, U>::Scalar>::value>::type * = nullptr>
0614     std::basic_ostream<char_t, traits_t> &operator<<(std::basic_ostream<char_t, traits_t> &os,
0615                                                      PositionVector3D<T, U> const &v)
0616     {
0617        if (os) {
0618           os << "{ ";
0619           for (std::size_t i = 0; i < PositionVector3D<T, U>::Scalar::size(); ++i) {
0620              os << "(" << v.x()[i] << "," << v.y()[i] << "," << v.z()[i] << ") ";
0621           }
0622           os << "}";
0623        }
0624        return os;
0625     } // op<< <>()
0626 
0627     template< class char_t, class traits_t, class T, class U >
0628       inline
0629       std::basic_istream<char_t,traits_t> &
0630       operator >> ( std::basic_istream<char_t,traits_t> & is
0631                   , PositionVector3D<T,U> & v
0632                   )
0633     {
0634       if( !is )  return is;
0635 
0636       typename T::Scalar a, b, c;
0637 
0638       if( detail::get_manip( is, detail::bitforbit ) )  {
0639         detail::set_manip( is, detail::bitforbit, '\00' );
0640         typedef GenVector_detail::BitReproducible BR;
0641         BR::Input(is, a);
0642         BR::Input(is, b);
0643         BR::Input(is, c);
0644       }
0645       else  {
0646         detail::require_delim( is, detail::open  );  is >> a;
0647         detail::require_delim( is, detail::sep   );  is >> b;
0648         detail::require_delim( is, detail::sep   );  is >> c;
0649         detail::require_delim( is, detail::close );
0650       }
0651 
0652       if( is )
0653         v.SetCoordinates(a, b, c);
0654       return is;
0655 
0656     }  // op>> <>()
0657 
0658 
0659 
0660 
0661   } // namespace Math
0662 
0663 } // namespace ROOT
0664 
0665 
0666 #endif /* ROOT_Math_GenVector_PositionVector3D  */