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     @sa Overview of the @ref GenVector "physics vector library"
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    PtEtaPhiE4D() : fPt(0), fEta(0), fPhi(0), fE(0) { }
0067
0068    /**
0069       Constructor  from pt, eta, phi, e values
0070    */
0071    PtEtaPhiE4D(Scalar pt, Scalar eta, Scalar phi, Scalar e) :
0072       fPt(pt), fEta(eta), fPhi(phi), fE(e) { Restrict(); }
0073
0074    /**
0075       Generic constructor from any 4D coordinate system implementing
0076       Pt(), Eta(), Phi() and E()
0077    */
0078    template <class CoordSystem >
0079    explicit constexpr PtEtaPhiE4D(const CoordSystem & c) :
0080       fPt(c.Pt()), fEta(c.Eta()), fPhi(c.Phi()), fE(c.E())  { }
0081
0082    // for g++  3.2 and 3.4 on 32 bits found that the compiler generated copy ctor and assignment are much slower
0083    // so we decided to re-implement them ( there is no no need to have them with g++4)
0084
0085    /**
0086       copy constructor
0087     */
0088    PtEtaPhiE4D(const PtEtaPhiE4D & v) :
0089       fPt(v.fPt), fEta(v.fEta), fPhi(v.fPhi), fE(v.fE) { }
0090
0091    /**
0092       assignment operator
0093     */
0094    PtEtaPhiE4D & operator = (const PtEtaPhiE4D & v) {
0095       fPt  = v.fPt;
0096       fEta = v.fEta;
0097       fPhi = v.fPhi;
0098       fE   = v.fE;
0099       return *this;
0100    }
0101
0102
0103    /**
0104       Set internal data based on an array of 4 Scalar numbers
0105    */
0106    void SetCoordinates( const Scalar src[] )
0107    { fPt=src[0]; fEta=src[1]; fPhi=src[2]; fE=src[3]; Restrict(); }
0108
0109    /**
0110       get internal data into an array of 4 Scalar numbers
0111    */
0112    void GetCoordinates( Scalar dest[] ) const
0113    { dest[0] = fPt; dest[1] = fEta; dest[2] = fPhi; dest[3] = fE; }
0114
0115    /**
0116       Set internal data based on 4 Scalar numbers
0117    */
0118    void SetCoordinates(Scalar pt, Scalar eta, Scalar phi, Scalar e)
0119    { fPt=pt; fEta = eta; fPhi = phi; fE = e; Restrict(); }
0120
0121    /**
0122       get internal data into 4 Scalar numbers
0123    */
0124    void
0125    GetCoordinates(Scalar& pt, Scalar & eta, Scalar & phi, Scalar& e) const
0126    { pt=fPt; eta=fEta; phi = fPhi; e = fE; }
0127
0128    // --------- Coordinates and Coordinate-like Scalar properties -------------
0129
0130    // 4-D Cylindrical eta coordinate accessors
0131
0132    Scalar Pt()  const { return fPt;  }
0133    Scalar Eta() const { return fEta; }
0134    Scalar Phi() const { return fPhi; }
0135    Scalar E()   const { return fE;   }
0136
0137    Scalar Perp()const { return Pt(); }
0138    Scalar Rho() const { return Pt(); }
0139    Scalar T()   const { return E();  }
0140
0141    // other coordinate representation
0142
0143    Scalar Px() const { using std::cos; return fPt * cos(fPhi); }
0144    Scalar X () const { return Px();         }
0145    Scalar Py() const { using std::sin; return fPt * sin(fPhi); }
0146    Scalar Y () const { return Py();         }
0147    Scalar Pz() const {
0148       using std:: sinh;
0149       return fPt > 0 ? fPt * sinh(fEta) : fEta == 0 ? 0 : fEta > 0 ? fEta - etaMax<Scalar>() : fEta + etaMax<Scalar>();
0150    }
0151    Scalar Z () const { return Pz(); }
0152
0153    /**
0154        magnitude of momentum
0155    */
0156    Scalar P() const {
0157      using std::cosh;
0158       return fPt > 0 ? fPt * cosh(fEta)
0159                      : fEta > etaMax<Scalar>() ? fEta - etaMax<Scalar>()
0160                                                : fEta < -etaMax<Scalar>() ? -fEta - etaMax<Scalar>() : 0;
0161    }
0162    Scalar R() const { return P(); }
0163
0164    /**
0165        squared magnitude of spatial components (momentum squared)
0166    */
0167    Scalar P2() const
0168    {
0169       const Scalar p = P();
0170       return p * p;
0171    }
0172
0173    /**
0174       vector magnitude squared (or mass squared)
0175    */
0176    Scalar M2() const
0177    {
0178       const Scalar p = P();
0179       return fE * fE - p * p;
0180    }
0181    Scalar Mag2() const { return M2(); }
0182
0183    /**
0184       invariant mass
0185    */
0186    Scalar M() const    {
0187       const Scalar mm = M2();
0188       if (mm >= 0) {
0189          using std::sqrt;
0190          return sqrt(mm);
0191       } else {
0192          GenVector::Throw ("PtEtaPhiE4D::M() - Tachyonic:\n"
0193                            "    Pt and Eta give P such that P^2 > E^2, so the mass would be imaginary");
0194          using std::sqrt;
0195          return -sqrt(-mm);
0196       }
0197    }
0198    Scalar Mag() const    { return M(); }
0199
0200    /**
0201        transverse spatial component squared
0202    */
0203    Scalar Pt2()   const { return fPt*fPt;}
0204    Scalar Perp2() const { return Pt2();  }
0205
0206    /**
0207        transverse mass squared
0208    */
0209    Scalar Mt2() const {  Scalar pz = Pz(); return fE*fE  - pz*pz; }
0210
0211    /**
0212       transverse mass
0213    */
0214    Scalar Mt() const {
0215       const Scalar mm = Mt2();
0216       if (mm >= 0) {
0217          using std::sqrt;
0218          return sqrt(mm);
0219       } else {
0220          GenVector::Throw ("PtEtaPhiE4D::Mt() - Tachyonic:\n"
0221                            "    Pt and Eta give Pz such that Pz^2 > E^2, so the mass would be imaginary");
0222          using std::sqrt;
0223          return -sqrt(-mm);
0224       }
0225    }
0226
0227    /**
0228       transverse energy
0229    */
0230    /**
0231       transverse energy
0232    */
0233    Scalar Et() const {
0234       using std::cosh;
0235       return fE / cosh(fEta); // faster using eta
0236    }
0237
0238    /**
0239        transverse energy squared
0240    */
0241    Scalar Et2() const
0242    {
0243       const Scalar et = Et();
0244       return et * et;
0245    }
0246
0247 private:
0248    inline static Scalar pi() { return M_PI; }
0249    inline void Restrict() {
0250       using std::floor;
0251       if (fPhi <= -pi() || fPhi > pi()) fPhi = fPhi - floor(fPhi / (2 * pi()) + .5) * 2 * pi();
0252    }
0253 public:
0254
0255    /**
0256       polar angle
0257    */
0258    Scalar Theta() const { using std::atan; return (fPt > 0 ? Scalar(2) * atan(exp(-fEta)) : fEta >= 0 ? 0 : pi()); }
0259
0260    // --------- Set Coordinates of this system  ---------------
0261
0262    /**
0263       set Pt value
0264    */
0265    void SetPt( Scalar  pt) {
0266       fPt = pt;
0267    }
0268    /**
0269       set eta value
0270    */
0271    void SetEta( Scalar  eta) {
0272       fEta = eta;
0273    }
0274    /**
0275       set phi value
0276    */
0277    void SetPhi( Scalar  phi) {
0278       fPhi = phi;
0279       Restrict();
0280    }
0281    /**
0282       set E value
0283    */
0284    void SetE( Scalar  e) {
0285       fE = e;
0286    }
0287
0288    /**
0289        set values using cartesian coordinate system
0290    */
0291    void SetPxPyPzE(Scalar px, Scalar py, Scalar pz, Scalar e);
0292
0293
0294    // ------ Manipulations -------------
0295
0296    /**
0297       negate the 4-vector
0298    */
0299    void Negate( ) {
0300       fPhi = ( fPhi > 0 ? fPhi - pi() : fPhi + pi()  );
0301       fEta = - fEta;
0302       fE = - fE;
0303    }
0304
0305    /**
0306       Scale coordinate values by a scalar quantity a
0307    */
0308    void Scale( Scalar a) {
0309       if (a < 0) {
0310          Negate(); a = -a;
0311       }
0312       fPt *= a;
0313       fE  *= a;
0314    }
0315
0316    /**
0317       Assignment from a generic coordinate system implementing
0318       Pt(), Eta(), Phi() and E()
0319    */
0320    template <class CoordSystem >
0321    PtEtaPhiE4D & operator = (const CoordSystem & c) {
0322       fPt  = c.Pt();
0323       fEta = c.Eta();
0324       fPhi = c.Phi();
0325       fE   = c.E();
0326       return *this;
0327    }
0328
0329    /**
0330       Exact equality
0331    */
0332    bool operator == (const PtEtaPhiE4D & rhs) const {
0333       return fPt == rhs.fPt && fEta == rhs.fEta
0334          && fPhi == rhs.fPhi && fE == rhs.fE;
0335    }
0336    bool operator != (const PtEtaPhiE4D & rhs) const {return !(operator==(rhs));}
0337
0338    // ============= Compatibility section ==================
0339
0340    // The following make this coordinate system look enough like a CLHEP
0341    // vector that an assignment member template can work with either
0342    Scalar x() const { return X(); }
0343    Scalar y() const { return Y(); }
0344    Scalar z() const { return Z(); }
0345    Scalar t() const { return E(); }
0346
0347
0348
0349 #if defined(__MAKECINT__) || defined(G__DICTIONARY)
0350
0351    // ====== Set member functions for coordinates in other systems =======
0352
0353    void SetPx(Scalar px);
0354
0355    void SetPy(Scalar py);
0356
0357    void SetPz(Scalar pz);
0358
0359    void SetM(Scalar m);
0360
0361
0362 #endif
0363
0364 private:
0365
0366    ScalarType fPt;
0367    ScalarType fEta;
0368    ScalarType fPhi;
0369    ScalarType fE;
0370
0371 };
0372
0373
0374 } // end namespace Math
0375 } // end namespace ROOT
0376
0377
0378
0379 // move implementations here to avoid circle dependencies
0380 #include "Math/GenVector/PxPyPzE4D.h"
0381 #if defined(__MAKECINT__) || defined(G__DICTIONARY)
0382 #include "Math/GenVector/PtEtaPhiM4D.h"
0383 #endif
0384
0385 namespace ROOT {
0386
0387 namespace Math {
0388
0389 template <class ScalarType>
0390 inline void PtEtaPhiE4D<ScalarType>::SetPxPyPzE(Scalar px, Scalar py, Scalar pz, Scalar e) {
0391    *this = PxPyPzE4D<Scalar> (px, py, pz, e);
0392 }
0393
0394
0395 #if defined(__MAKECINT__) || defined(G__DICTIONARY)
0396
0397   // ====== Set member functions for coordinates in other systems =======
0398
0399 template <class ScalarType>
0400 inline void PtEtaPhiE4D<ScalarType>::SetPx(Scalar px) {
0401    GenVector_exception e("PtEtaPhiE4D::SetPx() is not supposed to be called");
0402    throw e;
0403    PxPyPzE4D<Scalar> v(*this); v.SetPx(px); *this = PtEtaPhiE4D<Scalar>(v);
0404 }
0405 template <class ScalarType>
0406 inline void PtEtaPhiE4D<ScalarType>::SetPy(Scalar py) {
0407    GenVector_exception e("PtEtaPhiE4D::SetPx() is not supposed to be called");
0408    throw e;
0409    PxPyPzE4D<Scalar> v(*this); v.SetPy(py); *this = PtEtaPhiE4D<Scalar>(v);
0410 }
0411 template <class ScalarType>
0412 inline void PtEtaPhiE4D<ScalarType>::SetPz(Scalar pz) {
0413    GenVector_exception e("PtEtaPhiE4D::SetPx() is not supposed to be called");
0414    throw e;
0415    PxPyPzE4D<Scalar> v(*this); v.SetPz(pz); *this = PtEtaPhiE4D<Scalar>(v);
0416 }
0417 template <class ScalarType>
0418 inline void PtEtaPhiE4D<ScalarType>::SetM(Scalar m) {
0419    GenVector_exception e("PtEtaPhiE4D::SetM() is not supposed to be called");
0420    throw e;
0421    PtEtaPhiM4D<Scalar> v(*this); v.SetM(m);
0422    *this = PtEtaPhiE4D<Scalar>(v);
0423 }
0424
0425 #endif  // endif __MAKE__CINT || G__DICTIONARY
0426
0427 } // end namespace Math
0428
0429 } // end namespace ROOT
0430
0431
0432
0433
0434 #endif // ROOT_Math_GenVector_PtEtaPhiE4D