Math/GenVector/PtEtaPhiE4D.h

0001
0002 // @(#)root/mathcore:$Id$
0003 // Authors: W. Brown, M. Fischler, L. Moneta    2005
0004
0005 /**********************************************************************
0006 *                                                                    *
0007 * Copyright (c) 2005 , LCG ROOT FNAL MathLib Team                    *
0008 *                                                                    *
0009 *                                                                    *
0010 **********************************************************************/
0011
0012 // Header file for class PtEtaPhiE4D
0013 //
0014 // Created by: fischler at Wed Jul 20 2005
0015 //   based on CylindricalEta4D by moneta
0016 //
0017 // Last update: $Id$
0018 //
0019 #ifndef ROOT_Math_GenVector_PtEtaPhiE4D
0020 #define ROOT_Math_GenVector_PtEtaPhiE4D  1
0021
0022 #include "Math/Math.h"
0023
0024 #include "Math/GenVector/etaMax.h"
0025
0026 #include "Math/GenVector/GenVector_exception.h"
0027
0028
0029
0030 //#define TRACE_CE
0031 #ifdef TRACE_CE
0032 #include <iostream>
0033 #endif
0034
0035 #include <cmath>
0036
0037 namespace ROOT {
0038
0039 namespace Math {
0040
0041 //__________________________________________________________________________________________
0042 /**
0043     Class describing a 4D cylindrical coordinate system
0044     using Pt , Phi, Eta and E (or rho, phi, eta , T)
0045     The metric used is (-,-,-,+).
0046     Phi is restricted to be in the range [-PI,PI)
0047
0048     @ingroup GenVector
0049
0050     @see GenVector
0051 */
0052
0053 template <class ScalarType>
0054 class PtEtaPhiE4D {
0055
0056 public :
0057
0058    typedef ScalarType Scalar;
0059    static constexpr unsigned int Dimension = 4U;
0060
0061    // --------- Constructors ---------------
0062
0063    /**
0064       Default constructor gives zero 4-vector
0065    */
0066    constexpr PtEtaPhiE4D() noexcept = default;
0067
0068    /**
0069       Constructor  from pt, eta, phi, e values
0070    */
0071    constexpr PtEtaPhiE4D(Scalar pt, Scalar eta, Scalar phi, Scalar e) noexcept : fPt(pt), fEta(eta), fPhi(phi), fE(e)
0072    {
0073       Restrict();
0074    }
0075
0076    /**
0077       Generic constructor from any 4D coordinate system implementing
0078       Pt(), Eta(), Phi() and E()
0079    */
0080    template <class CoordSystem >
0081    explicit constexpr PtEtaPhiE4D(const CoordSystem & c) :
0082       fPt(c.Pt()), fEta(c.Eta()), fPhi(c.Phi()), fE(c.E())  { }
0083
0084    /**
0085       Set internal data based on an array of 4 Scalar numbers
0086    */
0087    void SetCoordinates( const Scalar src[] )
0088    { fPt=src[0]; fEta=src[1]; fPhi=src[2]; fE=src[3]; Restrict(); }
0089
0090    /**
0091       get internal data into an array of 4 Scalar numbers
0092    */
0093    void GetCoordinates( Scalar dest[] ) const
0094    { dest[0] = fPt; dest[1] = fEta; dest[2] = fPhi; dest[3] = fE; }
0095
0096    /**
0097       Set internal data based on 4 Scalar numbers
0098    */
0099    void SetCoordinates(Scalar pt, Scalar eta, Scalar phi, Scalar e)
0100    { fPt=pt; fEta = eta; fPhi = phi; fE = e; Restrict(); }
0101
0102    /**
0103       get internal data into 4 Scalar numbers
0104    */
0105    void
0106    GetCoordinates(Scalar& pt, Scalar & eta, Scalar & phi, Scalar& e) const
0107    { pt=fPt; eta=fEta; phi = fPhi; e = fE; }
0108
0109    // --------- Coordinates and Coordinate-like Scalar properties -------------
0110
0111    // 4-D Cylindrical eta coordinate accessors
0112
0113    Scalar Pt()  const { return fPt;  }
0114    Scalar Eta() const { return fEta; }
0115    Scalar Phi() const { return fPhi; }
0116    Scalar E()   const { return fE;   }
0117
0118    Scalar Perp()const { return Pt(); }
0119    Scalar Rho() const { return Pt(); }
0120    Scalar T()   const { return E();  }
0121
0122    // other coordinate representation
0123
0124    Scalar Px() const { using std::cos; return fPt * cos(fPhi); }
0125    Scalar X () const { return Px();         }
0126    Scalar Py() const { using std::sin; return fPt * sin(fPhi); }
0127    Scalar Y () const { return Py();         }
0128    Scalar Pz() const {
0129       using std:: sinh;
0130       return fPt > 0 ? fPt * sinh(fEta) : fEta == 0 ? 0 : fEta > 0 ? fEta - etaMax<Scalar>() : fEta + etaMax<Scalar>();
0131    }
0132    Scalar Z () const { return Pz(); }
0133
0134    /**
0135        magnitude of momentum
0136    */
0137    Scalar P() const {
0138      using std::cosh;
0139       return fPt > 0 ? fPt * cosh(fEta)
0140                      : fEta > etaMax<Scalar>() ? fEta - etaMax<Scalar>()
0141                                                : fEta < -etaMax<Scalar>() ? -fEta - etaMax<Scalar>() : 0;
0142    }
0143    Scalar R() const { return P(); }
0144
0145    /**
0146        squared magnitude of spatial components (momentum squared)
0147    */
0148    Scalar P2() const
0149    {
0150       const Scalar p = P();
0151       return p * p;
0152    }
0153
0154    /**
0155       vector magnitude squared (or mass squared)
0156    */
0157    Scalar M2() const
0158    {
0159       const Scalar p = P();
0160       return fE * fE - p * p;
0161    }
0162    Scalar Mag2() const { return M2(); }
0163
0164    /**
0165       invariant mass
0166    */
0167    Scalar M() const    {
0168       const Scalar mm = M2();
0169       if (mm >= 0) {
0170          using std::sqrt;
0171          return sqrt(mm);
0172       } else {
0173          GenVector_Throw("PtEtaPhiE4D::M() - Tachyonic:\n"
0174                          "    Pt and Eta give P such that P^2 > E^2, so the mass would be imaginary");
0175          using std::sqrt;
0176          return -sqrt(-mm);
0177       }
0178    }
0179    Scalar Mag() const    { return M(); }
0180
0181    /**
0182        transverse spatial component squared
0183    */
0184    Scalar Pt2()   const { return fPt*fPt;}
0185    Scalar Perp2() const { return Pt2();  }
0186
0187    /**
0188        transverse mass squared
0189    */
0190    Scalar Mt2() const {  Scalar pz = Pz(); return fE*fE  - pz*pz; }
0191
0192    /**
0193       transverse mass
0194    */
0195    Scalar Mt() const {
0196       const Scalar mm = Mt2();
0197       if (mm >= 0) {
0198          using std::sqrt;
0199          return sqrt(mm);
0200       } else {
0201          GenVector_Throw("PtEtaPhiE4D::Mt() - Tachyonic:\n"
0202                          "    Pt and Eta give Pz such that Pz^2 > E^2, so the mass would be imaginary");
0203          using std::sqrt;
0204          return -sqrt(-mm);
0205       }
0206    }
0207
0208    /**
0209       transverse energy
0210    */
0211    /**
0212       transverse energy
0213    */
0214    Scalar Et() const {
0215       using std::cosh;
0216       return fE / cosh(fEta); // faster using eta
0217    }
0218
0219    /**
0220        transverse energy squared
0221    */
0222    Scalar Et2() const
0223    {
0224       const Scalar et = Et();
0225       return et * et;
0226    }
0227
0228 private:
0229    inline static Scalar pi() { return M_PI; }
0230    inline void Restrict() {
0231       using std::floor;
0232       if (fPhi <= -pi() || fPhi > pi()) fPhi = fPhi - floor(fPhi / (2 * pi()) + .5) * 2 * pi();
0233    }
0234 public:
0235
0236    /**
0237       polar angle
0238    */
0239    Scalar Theta() const { using std::atan; return (fPt > 0 ? Scalar(2) * atan(exp(-fEta)) : fEta >= 0 ? 0 : pi()); }
0240
0241    // --------- Set Coordinates of this system  ---------------
0242
0243    /**
0244       set Pt value
0245    */
0246    void SetPt( Scalar  pt) {
0247       fPt = pt;
0248    }
0249    /**
0250       set eta value
0251    */
0252    void SetEta( Scalar  eta) {
0253       fEta = eta;
0254    }
0255    /**
0256       set phi value
0257    */
0258    void SetPhi( Scalar  phi) {
0259       fPhi = phi;
0260       Restrict();
0261    }
0262    /**
0263       set E value
0264    */
0265    void SetE( Scalar  e) {
0266       fE = e;
0267    }
0268
0269    /**
0270        set values using cartesian coordinate system
0271    */
0272    void SetPxPyPzE(Scalar px, Scalar py, Scalar pz, Scalar e);
0273
0274
0275    // ------ Manipulations -------------
0276
0277    /**
0278       negate the 4-vector
0279    */
0280    void Negate( ) {
0281       fPhi = ( fPhi > 0 ? fPhi - pi() : fPhi + pi()  );
0282       fEta = - fEta;
0283       fE = - fE;
0284    }
0285
0286    /**
0287       Scale coordinate values by a scalar quantity a
0288    */
0289    void Scale( Scalar a) {
0290       if (a < 0) {
0291          Negate(); a = -a;
0292       }
0293       fPt *= a;
0294       fE  *= a;
0295    }
0296
0297    /**
0298       Assignment from a generic coordinate system implementing
0299       Pt(), Eta(), Phi() and E()
0300    */
0301    template <class CoordSystem >
0302    PtEtaPhiE4D & operator = (const CoordSystem & c) {
0303       fPt  = c.Pt();
0304       fEta = c.Eta();
0305       fPhi = c.Phi();
0306       fE   = c.E();
0307       return *this;
0308    }
0309
0310    /**
0311       Exact equality
0312    */
0313    bool operator == (const PtEtaPhiE4D & rhs) const {
0314       return fPt == rhs.fPt && fEta == rhs.fEta
0315          && fPhi == rhs.fPhi && fE == rhs.fE;
0316    }
0317    bool operator != (const PtEtaPhiE4D & rhs) const {return !(operator==(rhs));}
0318
0319    // ============= Compatibility section ==================
0320
0321    // The following make this coordinate system look enough like a CLHEP
0322    // vector that an assignment member template can work with either
0323    Scalar x() const { return X(); }
0324    Scalar y() const { return Y(); }
0325    Scalar z() const { return Z(); }
0326    Scalar t() const { return E(); }
0327
0328
0329
0330 #if defined(__MAKECINT__) || defined(G__DICTIONARY)
0331
0332    // ====== Set member functions for coordinates in other systems =======
0333
0334    void SetPx(Scalar px);
0335
0336    void SetPy(Scalar py);
0337
0338    void SetPz(Scalar pz);
0339
0340    void SetM(Scalar m);
0341
0342 #endif
0343
0344 private:
0345    ScalarType fPt = 0.;
0346    ScalarType fEta = 0.;
0347    ScalarType fPhi = 0.;
0348    ScalarType fE = 0.;
0349 };
0350
0351
0352 } // end namespace Math
0353 } // end namespace ROOT
0354
0355
0356
0357 // move implementations here to avoid circle dependencies
0358 #include "Math/GenVector/PxPyPzE4D.h"
0359 #if defined(__MAKECINT__) || defined(G__DICTIONARY)
0360 #include "Math/GenVector/PtEtaPhiM4D.h"
0361 #endif
0362
0363 namespace ROOT {
0364
0365 namespace Math {
0366
0367 template <class ScalarType>
0368 inline void PtEtaPhiE4D<ScalarType>::SetPxPyPzE(Scalar px, Scalar py, Scalar pz, Scalar e) {
0369    *this = PxPyPzE4D<Scalar> (px, py, pz, e);
0370 }
0371
0372
0373 #if defined(__MAKECINT__) || defined(G__DICTIONARY)
0374
0375   // ====== Set member functions for coordinates in other systems =======
0376
0377 template <class ScalarType>
0378 inline void PtEtaPhiE4D<ScalarType>::SetPx(Scalar px) {
0379    GenVector_exception e("PtEtaPhiE4D::SetPx() is not supposed to be called");
0380    throw e;
0381    PxPyPzE4D<Scalar> v(*this); v.SetPx(px); *this = PtEtaPhiE4D<Scalar>(v);
0382 }
0383 template <class ScalarType>
0384 inline void PtEtaPhiE4D<ScalarType>::SetPy(Scalar py) {
0385    GenVector_exception e("PtEtaPhiE4D::SetPx() is not supposed to be called");
0386    throw e;
0387    PxPyPzE4D<Scalar> v(*this); v.SetPy(py); *this = PtEtaPhiE4D<Scalar>(v);
0388 }
0389 template <class ScalarType>
0390 inline void PtEtaPhiE4D<ScalarType>::SetPz(Scalar pz) {
0391    GenVector_exception e("PtEtaPhiE4D::SetPx() is not supposed to be called");
0392    throw e;
0393    PxPyPzE4D<Scalar> v(*this); v.SetPz(pz); *this = PtEtaPhiE4D<Scalar>(v);
0394 }
0395 template <class ScalarType>
0396 inline void PtEtaPhiE4D<ScalarType>::SetM(Scalar m) {
0397    GenVector_exception e("PtEtaPhiE4D::SetM() is not supposed to be called");
0398    throw e;
0399    PtEtaPhiM4D<Scalar> v(*this); v.SetM(m);
0400    *this = PtEtaPhiE4D<Scalar>(v);
0401 }
0402
0403 #endif  // endif __MAKE__CINT || G__DICTIONARY
0404
0405 } // end namespace Math
0406
0407 } // end namespace ROOT
0408
0409
0410
0411
0412 #endif // ROOT_Math_GenVector_PtEtaPhiE4D