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 //   Successor 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 #include "TMath.h"
0023
0024 #include "Math/GenVector/etaMax.h"
0025
0026 #include "Math/GenVector/GenVector_exception.h"
0027
0028
0029 //#define TRACE_CE
0030 #ifdef TRACE_CE
0031 #include <iostream>
0032 #endif
0033
0034 #include <cmath>
0035
0036 namespace ROOT {
0037
0038 namespace Math {
0039
0040 //__________________________________________________________________________________________
0041 /**
0042     Class describing a 4D cylindrical coordinate system
0043     using Pt , Phi, Eta and M (mass)
0044     The metric used is (-,-,-,+).
0045     Spacelike particles (M2 < 0) are described with negative mass values,
0046     but in this case m2 must always be less than P2 to preserve a positive value of E2
0047     Phi is restricted to be in the range [-PI,PI)
0048
0049     @ingroup GenVector
0050
0051     @see GenVector
0052 */
0053
0054 template <class ScalarType>
0055 class PtEtaPhiM4D {
0056
0057 public :
0058
0059    typedef ScalarType Scalar;
0060    static constexpr unsigned int Dimension = 4U;
0061
0062    // --------- Constructors ---------------
0063
0064    /**
0065       Default constructor gives zero 4-vector (with zero mass)
0066    */
0067    PtEtaPhiM4D() : fPt(0), fEta(0), fPhi(0), fM(0) { }
0068
0069    /**
0070       Constructor  from pt, eta, phi, mass values
0071    */
0072    PtEtaPhiM4D(Scalar pt, Scalar eta, Scalar phi, Scalar mass) :
0073       fPt(pt), fEta(eta), fPhi(phi), fM(mass) {
0074       RestrictPhi();
0075       if (fM < 0) RestrictNegMass();
0076    }
0077
0078    /**
0079       Generic constructor from any 4D coordinate system implementing
0080       Pt(), Eta(), Phi() and M()
0081    */
0082    template <class CoordSystem >
0083    explicit constexpr PtEtaPhiM4D(const CoordSystem & c) :
0084       fPt(c.Pt()), fEta(c.Eta()), fPhi(c.Phi()), fM(c.M())  { RestrictPhi(); }
0085
0086
0087    /**
0088       Set internal data based on an array of 4 Scalar numbers
0089    */
0090    void SetCoordinates( const Scalar src[] ) {
0091       fPt=src[0]; fEta=src[1]; fPhi=src[2]; fM=src[3];
0092       RestrictPhi();
0093       if (fM <0) RestrictNegMass();
0094    }
0095
0096    /**
0097       get internal data into an array of 4 Scalar numbers
0098    */
0099    void GetCoordinates( Scalar dest[] ) const
0100    { dest[0] = fPt; dest[1] = fEta; dest[2] = fPhi; dest[3] = fM; }
0101
0102    /**
0103       Set internal data based on 4 Scalar numbers
0104    */
0105    void SetCoordinates(Scalar pt, Scalar eta, Scalar phi, Scalar mass) {
0106       fPt=pt; fEta = eta; fPhi = phi; fM = mass;
0107       RestrictPhi();
0108       if (fM <0) RestrictNegMass();
0109    }
0110
0111    /**
0112       get internal data into 4 Scalar numbers
0113    */
0114    void
0115    GetCoordinates(Scalar& pt, Scalar & eta, Scalar & phi, Scalar& mass) const
0116    { pt=fPt; eta=fEta; phi = fPhi; mass = fM; }
0117
0118    // --------- Coordinates and Coordinate-like Scalar properties -------------
0119
0120    // 4-D Cylindrical eta coordinate accessors
0121
0122    Scalar Pt()  const { return fPt;  }
0123    Scalar Eta() const { return fEta; }
0124    Scalar Phi() const { return fPhi; }
0125    /**
0126        M() is the invariant mass;
0127        in this coordinate system it can be negagative if set that way.
0128    */
0129    Scalar M()   const { return fM;   }
0130    Scalar Mag() const { return M(); }
0131
0132    Scalar Perp()const { return Pt(); }
0133    Scalar Rho() const { return Pt(); }
0134
0135    // other coordinate representation
0136
0137    Scalar Px() const { using std::cos; return fPt * cos(fPhi); }
0138    Scalar X () const { return Px();         }
0139    Scalar Py() const { using std::sin; return fPt * sin(fPhi); }
0140    Scalar Y () const { return Py();         }
0141    Scalar Pz() const {
0142       using std::sinh;
0143       return fPt > 0 ? fPt * sinh(fEta) : fEta == 0 ? 0 : fEta > 0 ? fEta - etaMax<Scalar>() : fEta + etaMax<Scalar>();
0144    }
0145    Scalar Z () const { return Pz(); }
0146
0147    /**
0148        magnitude of momentum
0149    */
0150    Scalar P() const {
0151       using std::cosh;
0152       return fPt > 0 ? fPt * cosh(fEta)
0153                      : fEta > etaMax<Scalar>() ? fEta - etaMax<Scalar>()
0154                                                : fEta < -etaMax<Scalar>() ? -fEta - etaMax<Scalar>() : 0;
0155    }
0156    Scalar R() const { return P(); }
0157
0158    /**
0159        squared magnitude of spatial components (momentum squared)
0160    */
0161    Scalar P2() const
0162    {
0163       const Scalar p = P();
0164       return p * p;
0165    }
0166
0167    /**
0168        energy squared
0169    */
0170    Scalar E2() const {
0171       Scalar e2 =  P2() + M2();
0172       // avoid rounding error which can make E2 negative when M2 is negative
0173       return e2 > 0 ? e2 : 0;
0174    }
0175
0176    /**
0177        Energy (timelike component of momentum-energy 4-vector)
0178    */
0179    Scalar E() const { using std::sqrt; return sqrt(E2()); }
0180
0181    Scalar T()   const { return E();  }
0182
0183    /**
0184       vector magnitude squared (or mass squared)
0185       In case of negative mass (spacelike particles return negative values)
0186    */
0187    Scalar M2() const   {
0188       return ( fM  >= 0 ) ?  fM*fM :  -fM*fM;
0189    }
0190    Scalar Mag2() const { return M2();  }
0191
0192    /**
0193        transverse spatial component squared
0194    */
0195    Scalar Pt2()   const { return fPt*fPt;}
0196    Scalar Perp2() const { return Pt2();  }
0197
0198    /**
0199        transverse mass squared
0200    */
0201    Scalar Mt2() const { return M2()  + fPt*fPt; }
0202
0203    /**
0204       transverse mass - will be negative if Mt2() is negative
0205    */
0206    Scalar Mt() const {
0207       const Scalar mm = Mt2();
0208       if (mm >= 0) {
0209          using std::sqrt;
0210          return sqrt(mm);
0211       } else {
0212          GenVector_Throw("PtEtaPhiM4D::Mt() - Tachyonic:\n"
0213                          "    Pz^2 > E^2 so the transverse mass would be imaginary");
0214          using std::sqrt;
0215          return -sqrt(-mm);
0216       }
0217    }
0218
0219    /**
0220        transverse energy squared
0221    */
0222    Scalar Et2() const {
0223       // a bit faster than et * et
0224       using std::cosh;
0225       return 2. * E2() / (cosh(2 * fEta) + 1);
0226    }
0227
0228    /**
0229       transverse energy
0230    */
0231    Scalar Et() const { using std::cosh; return E() / cosh(fEta); }
0232
0233 private:
0234    inline static Scalar pi() { return TMath::Pi(); }
0235    inline void RestrictPhi() {
0236       using std::floor;
0237       if (fPhi <= -pi() || fPhi > pi()) fPhi = fPhi - floor(fPhi / (2 * pi()) + .5) * 2 * pi();
0238    }
0239    // restrict the value of negative mass to avoid unphysical negative E2 values
0240    // M2 must be less than P2 for the tachionic particles - otherwise use positive values
0241    inline void RestrictNegMass() {
0242       if (fM < 0) {
0243          if (P2() - fM * fM < 0) {
0244             GenVector_Throw("PtEtaPhiM4D::unphysical value of mass, set to closest physical value");
0245             fM = -P();
0246          }
0247       }
0248    }
0249
0250 public:
0251
0252    /**
0253       polar angle
0254    */
0255    Scalar Theta() const { using std::atan; return (fPt > 0 ? Scalar(2) * atan(exp(-fEta)) : fEta >= 0 ? 0 : pi()); }
0256
0257    // --------- Set Coordinates of this system  ---------------
0258
0259    /**
0260       set Pt value
0261    */
0262    void SetPt( Scalar  pt) {
0263       fPt = pt;
0264    }
0265    /**
0266       set eta value
0267    */
0268    void SetEta( Scalar  eta) {
0269       fEta = eta;
0270    }
0271    /**
0272       set phi value
0273    */
0274    void SetPhi( Scalar  phi) {
0275       fPhi = phi;
0276       RestrictPhi();
0277    }
0278    /**
0279       set M value
0280    */
0281    void SetM( Scalar  mass) {
0282       fM = mass;
0283       if (fM <0) RestrictNegMass();
0284    }
0285
0286    /**
0287        set values using cartesian coordinate system
0288    */
0289    void SetPxPyPzE(Scalar px, Scalar py, Scalar pz, Scalar e);
0290
0291
0292    // ------ Manipulations -------------
0293
0294    /**
0295       negate the 4-vector -- Note that the energy cannot be negate (would need an additional data member)
0296       therefore negate will work only on the spatial components
0297       One would need to use negate only with vectors having the energy as data members
0298    */
0299    void Negate( ) {
0300       fPhi = ( (fPhi > 0) ? fPhi - pi() : fPhi + pi()  );
0301       fEta = - fEta;
0302       GenVector_Throw("PtEtaPhiM4D::Negate - cannot negate the energy - can negate only the spatial components");
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       fM  *= a;
0314    }
0315
0316    /**
0317       Assignment from a generic coordinate system implementing
0318       Pt(), Eta(), Phi() and M()
0319    */
0320    template <class CoordSystem >
0321    PtEtaPhiM4D & operator = (const CoordSystem & c) {
0322       fPt  = c.Pt();
0323       fEta = c.Eta();
0324       fPhi = c.Phi();
0325       fM   = c.M();
0326       return *this;
0327    }
0328
0329    /**
0330       Exact equality
0331    */
0332    bool operator == (const PtEtaPhiM4D & rhs) const {
0333       return fPt == rhs.fPt && fEta == rhs.fEta
0334          && fPhi == rhs.fPhi && fM == rhs.fM;
0335    }
0336    bool operator != (const PtEtaPhiM4D & 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 #if defined(__MAKECINT__) || defined(G__DICTIONARY)
0349
0350    // ====== Set member functions for coordinates in other systems =======
0351
0352    void SetPx(Scalar px);
0353
0354    void SetPy(Scalar py);
0355
0356    void SetPz(Scalar pz);
0357
0358    void SetE(Scalar t);
0359
0360 #endif
0361
0362 private:
0363
0364    ScalarType fPt;
0365    ScalarType fEta;
0366    ScalarType fPhi;
0367    ScalarType fM;
0368
0369 };
0370
0371
0372 } // end namespace Math
0373 } // end namespace ROOT
0374
0375
0376 // move implementations here to avoid circle dependencies
0377 #include "Math/GenVector/PxPyPzE4D.h"
0378
0379
0380
0381 namespace ROOT {
0382
0383 namespace Math {
0384
0385
0386 template <class ScalarType>
0387 inline void PtEtaPhiM4D<ScalarType>::SetPxPyPzE(Scalar px, Scalar py, Scalar pz, Scalar e) {
0388    *this = PxPyPzE4D<Scalar> (px, py, pz, e);
0389 }
0390
0391
0392 #if defined(__MAKECINT__) || defined(G__DICTIONARY)
0393
0394   // ====== Set member functions for coordinates in other systems =======
0395
0396 template <class ScalarType>
0397 void PtEtaPhiM4D<ScalarType>::SetPx(Scalar px) {
0398    GenVector_exception e("PtEtaPhiM4D::SetPx() is not supposed to be called");
0399    throw e;
0400    PxPyPzE4D<Scalar> v(*this); v.SetPx(px); *this = PtEtaPhiM4D<Scalar>(v);
0401 }
0402 template <class ScalarType>
0403 void PtEtaPhiM4D<ScalarType>::SetPy(Scalar py) {
0404    GenVector_exception e("PtEtaPhiM4D::SetPx() is not supposed to be called");
0405    throw e;
0406    PxPyPzE4D<Scalar> v(*this); v.SetPy(py); *this = PtEtaPhiM4D<Scalar>(v);
0407 }
0408 template <class ScalarType>
0409 void PtEtaPhiM4D<ScalarType>::SetPz(Scalar pz) {
0410    GenVector_exception e("PtEtaPhiM4D::SetPx() is not supposed to be called");
0411    throw e;
0412    PxPyPzE4D<Scalar> v(*this); v.SetPz(pz); *this = PtEtaPhiM4D<Scalar>(v);
0413 }
0414 template <class ScalarType>
0415 void PtEtaPhiM4D<ScalarType>::SetE(Scalar energy) {
0416    GenVector_exception e("PtEtaPhiM4D::SetE() is not supposed to be called");
0417    throw e;
0418    PxPyPzE4D<Scalar> v(*this); v.SetE(energy);   *this = PtEtaPhiM4D<Scalar>(v);
0419 }
0420
0421 #endif  // endif __MAKE__CINT || G__DICTIONARY
0422
0423 } // end namespace Math
0424
0425 } // end namespace ROOT
0426
0427
0428
0429 #endif // ROOT_Math_GenVector_PtEtaPhiM4D