Math/GenVector/PtEtaPhiM4D.h

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