Math/GenVector/CylindricalEta3D.h

0001 // @(#)root/mathcore:$Id$
0002 // Authors: W. Brown, M. Fischler, L. Moneta    2005
0003
0004  /**********************************************************************
0005   *                                                                    *
0006   * Copyright (c) 2005 , LCG ROOT MathLib Team  and                    *
0007   *                      FNAL LCG ROOT MathLib Team                    *
0008   *                                                                    *
0009   *                                                                    *
0010   **********************************************************************/
0011
0012 // Header file for class CylindricalEta3D
0013 //
0014 // Created by: Lorenzo Moneta  at Mon May 30 11:58:46 2005
0015 // Major revamp:  M. Fischler  at Fri Jun 10 2005
0016 //
0017 // Last update: $Id$
0018
0019 //
0020 #ifndef ROOT_Math_GenVector_CylindricalEta3D
0021 #define ROOT_Math_GenVector_CylindricalEta3D  1
0022
0023 #include "Math/Math.h"
0024
0025 #include "Math/GenVector/etaMax.h"
0026
0027
0028 #include <limits>
0029 #include <cmath>
0030
0031
0032 namespace ROOT {
0033
0034 namespace Math {
0035
0036 //__________________________________________________________________________________________
0037   /**
0038       Class describing a cylindrical coordinate system based on eta (pseudorapidity) instead of z.
0039       The base coordinates are rho (transverse component) , eta and phi
0040       Phi is restricted to be in the range [-PI,PI)
0041
0042       @ingroup GenVector
0043
0044       @sa Overview of the @ref GenVector "physics vector library"
0045   */
0046
0047 template <class T>
0048 class CylindricalEta3D {
0049
0050 public :
0051
0052   typedef T Scalar;
0053   static constexpr unsigned int Dimension = 3U;
0054
0055   /**
0056      Default constructor with rho=eta=phi=0
0057    */
0058   CylindricalEta3D() : fRho(0), fEta(0), fPhi(0) {  }
0059
0060   /**
0061      Construct from rho eta and phi values
0062    */
0063   CylindricalEta3D(Scalar rho, Scalar eta, Scalar phi) :
0064              fRho(rho), fEta(eta), fPhi(phi) { Restrict(); }
0065
0066   /**
0067      Construct from any Vector or coordinate system implementing
0068      Rho(), Eta() and Phi()
0069     */
0070   template <class CoordSystem >
0071   explicit CylindricalEta3D( const CoordSystem & v ) :
0072      fRho(v.Rho() ),  fEta(v.Eta() ),  fPhi(v.Phi() )
0073   {
0074      using std::log;
0075      static Scalar bigEta = Scalar(-0.3) * log(std::numeric_limits<Scalar>::epsilon());
0076      if (std::fabs(fEta) > bigEta) {
0077         // This gives a small absolute adjustment in rho,
0078         // which, for large eta, results in a significant
0079         // improvement in the faithfullness of reproducing z.
0080         fRho *= v.Z() / Z();
0081     }
0082   }
0083
0084    // for g++  3.2 and 3.4 on 32 bits found that the compiler generated copy ctor and assignment are much slower
0085    // re-implement them ( there is no no need to have them with g++4)
0086
0087    /**
0088       copy constructor
0089    */
0090    CylindricalEta3D(const CylindricalEta3D & v) :
0091       fRho(v.Rho() ),  fEta(v.Eta() ),  fPhi(v.Phi() )  {   }
0092
0093    /**
0094       assignment operator
0095    */
0096    CylindricalEta3D & operator= (const CylindricalEta3D & v) {
0097       fRho = v.Rho();
0098       fEta = v.Eta();
0099       fPhi = v.Phi();
0100       return *this;
0101    }
0102
0103    /**
0104       Set internal data based on an array of 3 Scalar numbers
0105    */
0106    void SetCoordinates( const Scalar src[] )
0107    { fRho=src[0]; fEta=src[1]; fPhi=src[2]; Restrict(); }
0108
0109    /**
0110       get internal data into an array of 3 Scalar numbers
0111    */
0112    void GetCoordinates( Scalar dest[] ) const
0113    { dest[0] = fRho; dest[1] = fEta; dest[2] = fPhi; }
0114
0115    /**
0116       Set internal data based on 3 Scalar numbers
0117    */
0118    void SetCoordinates(Scalar  rho, Scalar  eta, Scalar  phi)
0119    { fRho=rho; fEta=eta; fPhi=phi; Restrict(); }
0120
0121    /**
0122       get internal data into 3 Scalar numbers
0123    */
0124    void GetCoordinates(Scalar& rho, Scalar& eta, Scalar& phi) const
0125    {rho=fRho; eta=fEta; phi=fPhi;}
0126
0127 private:
0128    inline static Scalar pi() { return M_PI; }
0129    inline void Restrict() {
0130       using std::floor;
0131       if (fPhi <= -pi() || fPhi > pi()) fPhi = fPhi - floor(fPhi / (2 * pi()) + .5) * 2 * pi();
0132       return;
0133    }
0134 public:
0135
0136    // accessors
0137
0138    T Rho()   const { return fRho; }
0139    T Eta()   const { return fEta; }
0140    T Phi()   const { return fPhi; }
0141    T X() const { using std::cos; return fRho * cos(fPhi); }
0142    T Y() const { using std::sin; return fRho * sin(fPhi); }
0143    T Z() const
0144    {
0145       using std::sinh;
0146       return fRho > 0 ? fRho * sinh(fEta) : fEta == 0 ? 0 : fEta > 0 ? fEta - etaMax<T>() : fEta + etaMax<T>();
0147    }
0148    T R() const
0149    {
0150       using std::cosh;
0151       return fRho > 0 ? fRho * cosh(fEta)
0152                       : fEta > etaMax<T>() ? fEta - etaMax<T>() : fEta < -etaMax<T>() ? -fEta - etaMax<T>() : 0;
0153    }
0154    T Mag2() const
0155    {
0156       const Scalar r = R();
0157       return r * r;
0158    }
0159    T Perp2() const { return fRho*fRho;            }
0160    T Theta() const { using std::atan; return fRho > 0 ? 2 * atan(exp(-fEta)) : (fEta >= 0 ? 0 : pi()); }
0161
0162    // setters (only for data members)
0163
0164
0165    /**
0166        set the rho coordinate value keeping eta and phi constant
0167    */
0168    void SetRho(T rho) {
0169       fRho = rho;
0170    }
0171
0172    /**
0173        set the eta coordinate value keeping rho and phi constant
0174    */
0175    void SetEta(T eta) {
0176       fEta = eta;
0177    }
0178
0179    /**
0180        set the phi coordinate value keeping rho and eta constant
0181    */
0182    void SetPhi(T phi) {
0183       fPhi = phi;
0184       Restrict();
0185    }
0186
0187    /**
0188        set all values using cartesian coordinates
0189    */
0190    void SetXYZ(Scalar x, Scalar y, Scalar z);
0191
0192
0193    /**
0194       scale by a scalar quantity a --
0195       for cylindrical eta coords, as long as a >= 0, only rho changes!
0196    */
0197    void Scale (T a) {
0198       if (a < 0) {
0199          Negate();
0200          a = -a;
0201       }
0202       // angles do not change when scaling by a positive quantity
0203       if (fRho > 0) {
0204          fRho *= a;
0205       } else if ( fEta >  etaMax<T>() ) {
0206          fEta =  ( fEta-etaMax<T>())*a + etaMax<T>();
0207       } else if ( fEta < -etaMax<T>() ) {
0208          fEta =  ( fEta+etaMax<T>())*a - etaMax<T>();
0209       } // when rho==0 and eta is not above etaMax, vector represents 0
0210       // and remains unchanged
0211    }
0212
0213    /**
0214       negate the vector
0215    */
0216    void Negate ( ) {
0217       fPhi = ( fPhi > 0 ? fPhi - pi() : fPhi + pi()  );
0218       fEta = -fEta;
0219    }
0220
0221    // assignment operators
0222    /**
0223       generic assignment operator from any coordinate system
0224    */
0225    template <class CoordSystem >
0226    CylindricalEta3D & operator= ( const CoordSystem & c ) {
0227       fRho = c.Rho();
0228       fEta = c.Eta();
0229       fPhi = c.Phi();
0230       return *this;
0231    }
0232
0233    /**
0234       Exact component-by-component equality
0235       Note: Peculiar representations of the zero vector such as (0,1,0) will
0236       not test as equal to one another.
0237    */
0238    bool operator==(const CylindricalEta3D & rhs) const {
0239       return fRho == rhs.fRho && fEta == rhs.fEta && fPhi == rhs.fPhi;
0240    }
0241    bool operator!= (const CylindricalEta3D & rhs) const
0242    {return !(operator==(rhs));}
0243
0244
0245    // ============= Compatibility section ==================
0246
0247    // The following make this coordinate system look enough like a CLHEP
0248    // vector that an assignment member template can work with either
0249    T x() const { return X();}
0250    T y() const { return Y();}
0251    T z() const { return Z(); }
0252
0253    // ============= Specializations for improved speed ==================
0254
0255    // (none)
0256
0257 #if defined(__MAKECINT__) || defined(G__DICTIONARY)
0258
0259    // ====== Set member functions for coordinates in other systems =======
0260
0261    void SetX(Scalar x);
0262
0263    void SetY(Scalar y);
0264
0265    void SetZ(Scalar z);
0266
0267    void SetR(Scalar r);
0268
0269    void SetTheta(Scalar theta);
0270
0271
0272 #endif
0273
0274
0275 private:
0276    T fRho;
0277    T fEta;
0278    T fPhi;
0279
0280 };
0281
0282   } // end namespace Math
0283
0284 } // end namespace ROOT
0285
0286
0287 // move implementations here to avoid circle dependencies
0288
0289 #include "Math/GenVector/Cartesian3D.h"
0290
0291 #if defined(__MAKECINT__) || defined(G__DICTIONARY)
0292 #include "Math/GenVector/GenVector_exception.h"
0293 #include "Math/GenVector/Polar3D.h"
0294 #endif
0295
0296 namespace ROOT {
0297
0298   namespace Math {
0299
0300 template <class T>
0301 void CylindricalEta3D<T>::SetXYZ(Scalar xx, Scalar yy, Scalar zz) {
0302    *this = Cartesian3D<Scalar>(xx, yy, zz);
0303 }
0304
0305 #if defined(__MAKECINT__) || defined(G__DICTIONARY)
0306
0307
0308      // ====== Set member functions for coordinates in other systems =======
0309
0310
0311 template <class T>
0312 void CylindricalEta3D<T>::SetX(Scalar xx) {
0313    GenVector_exception e("CylindricalEta3D::SetX() is not supposed to be called");
0314    throw e;
0315    Cartesian3D<Scalar> v(*this); v.SetX(xx);
0316    *this = CylindricalEta3D<Scalar>(v);
0317 }
0318 template <class T>
0319 void CylindricalEta3D<T>::SetY(Scalar yy) {
0320    GenVector_exception e("CylindricalEta3D::SetY() is not supposed to be called");
0321    throw e;
0322    Cartesian3D<Scalar> v(*this); v.SetY(yy);
0323    *this = CylindricalEta3D<Scalar>(v);
0324 }
0325 template <class T>
0326 void CylindricalEta3D<T>::SetZ(Scalar zz) {
0327    GenVector_exception e("CylindricalEta3D::SetZ() is not supposed to be called");
0328    throw e;
0329    Cartesian3D<Scalar> v(*this); v.SetZ(zz);
0330    *this = CylindricalEta3D<Scalar>(v);
0331 }
0332 template <class T>
0333 void CylindricalEta3D<T>::SetR(Scalar r) {
0334    GenVector_exception e("CylindricalEta3D::SetR() is not supposed to be called");
0335    throw e;
0336    Polar3D<Scalar> v(*this); v.SetR(r);
0337    *this = CylindricalEta3D<Scalar>(v);
0338 }
0339 template <class T>
0340 void CylindricalEta3D<T>::SetTheta(Scalar theta) {
0341    GenVector_exception e("CylindricalEta3D::SetTheta() is not supposed to be called");
0342    throw e;
0343    Polar3D<Scalar> v(*this); v.SetTheta(theta);
0344    *this = CylindricalEta3D<Scalar>(v);
0345 }
0346
0347 #endif
0348
0349
0350   } // end namespace Math
0351
0352 } // end namespace ROOT
0353
0354
0355
0356 #endif /* ROOT_Math_GenVector_CylindricalEta3D  */