// @(#)root/mathcore:$Id: PtEtaPhiM4D.h 24923 2008-07-23 15:43:05Z moneta $ // Authors: W. Brown, M. Fischler, L. Moneta 2005 /********************************************************************** * * * Copyright (c) 2005 , LCG ROOT FNAL MathLib Team * * * * * **********************************************************************/ // Header file for class PtEtaPhiM4D // // Created by: fischler at Wed Jul 21 2005 // Similar to PtEtaPhiMSystem by moneta // // Last update: $Id: PtEtaPhiM4D.h 24923 2008-07-23 15:43:05Z moneta $ // #ifndef ROOT_Math_GenVector_PtEtaPhiM4D #define ROOT_Math_GenVector_PtEtaPhiM4D 1 #ifndef ROOT_Math_Math #include "Math/Math.h" #endif #ifndef ROOT_Math_GenVector_etaMax #include "Math/GenVector/etaMax.h" #endif #ifndef ROOT_Math_GenVector_GenVector_exception #include "Math/GenVector/GenVector_exception.h" #endif //#define TRACE_CE #ifdef TRACE_CE #include #endif namespace ROOT { namespace Math { //__________________________________________________________________________________________ /** Class describing a 4D cylindrical coordinate system using Pt , Phi, Eta and M (mass) The metric used is (-,-,-,+). Spacelike particles (M2 < 0) are described with negative mass values, but in this case m2 must alwasy be less than P2 to preserve a positive value of E2 Phi is restricted to be in the range [-PI,PI) @ingroup GenVector */ template class PtEtaPhiM4D { public : typedef ScalarType Scalar; // --------- Constructors --------------- /** Default constructor gives zero 4-vector (with zero mass) */ PtEtaPhiM4D() : fPt(0), fEta(0), fPhi(0), fM(0) { } /** Constructor from pt, eta, phi, mass values */ PtEtaPhiM4D(Scalar pt, Scalar eta, Scalar phi, Scalar mass) : fPt(pt), fEta(eta), fPhi(phi), fM(mass) { RestrictPhi(); if (fM < 0) RestrictNegMass(); } /** Generic constructor from any 4D coordinate system implementing Pt(), Eta(), Phi() and M() */ template explicit PtEtaPhiM4D(const CoordSystem & c) : fPt(c.Pt()), fEta(c.Eta()), fPhi(c.Phi()), fM(c.M()) { RestrictPhi(); } // for g++ 3.2 and 3.4 on 32 bits found that the compiler generated copy ctor and assignment are much slower // so we decided to re-implement them ( there is no no need to have them with g++4) /** copy constructor */ PtEtaPhiM4D(const PtEtaPhiM4D & v) : fPt(v.fPt), fEta(v.fEta), fPhi(v.fPhi), fM(v.fM) { } /** assignment operator */ PtEtaPhiM4D & operator = (const PtEtaPhiM4D & v) { fPt = v.fPt; fEta = v.fEta; fPhi = v.fPhi; fM = v.fM; return *this; } /** Set internal data based on an array of 4 Scalar numbers */ void SetCoordinates( const Scalar src[] ) { fPt=src[0]; fEta=src[1]; fPhi=src[2]; fM=src[3]; RestrictPhi(); if (fM <0) RestrictNegMass(); } /** get internal data into an array of 4 Scalar numbers */ void GetCoordinates( Scalar dest[] ) const { dest[0] = fPt; dest[1] = fEta; dest[2] = fPhi; dest[3] = fM; } /** Set internal data based on 4 Scalar numbers */ void SetCoordinates(Scalar pt, Scalar eta, Scalar phi, Scalar mass) { fPt=pt; fEta = eta; fPhi = phi; fM = mass; RestrictPhi(); if (fM <0) RestrictNegMass(); } /** get internal data into 4 Scalar numbers */ void GetCoordinates(Scalar& pt, Scalar & eta, Scalar & phi, Scalar& mass) const { pt=fPt; eta=fEta; phi = fPhi; mass = fM; } // --------- Coordinates and Coordinate-like Scalar properties ------------- // 4-D Cylindrical eta coordinate accessors Scalar Pt() const { return fPt; } Scalar Eta() const { return fEta; } Scalar Phi() const { return fPhi; } /** M() is the invariant mass; in this coordinate system it can be negagative if set that way. */ Scalar M() const { return fM; } Scalar Mag() const { return M(); } Scalar Perp()const { return Pt(); } Scalar Rho() const { return Pt(); } // other coordinate representation Scalar Px() const { return fPt*cos(fPhi);} Scalar X () const { return Px(); } Scalar Py() const { return fPt*sin(fPhi);} Scalar Y () const { return Py(); } Scalar Pz() const { return fPt > 0 ? fPt*std::sinh(fEta) : fEta == 0 ? 0 : fEta > 0 ? fEta - etaMax() : fEta + etaMax() ; } Scalar Z () const { return Pz(); } /** magnitude of momentum */ Scalar P() const { return fPt > 0 ? fPt*std::cosh(fEta) : fEta > etaMax() ? fEta - etaMax() : fEta < -etaMax() ? -fEta - etaMax() : 0 ; } Scalar R() const { return P(); } /** squared magnitude of spatial components (momentum squared) */ Scalar P2() const { Scalar p = P(); return p*p; } /** energy squared */ Scalar E2() const { Scalar e2 = P2() + M2(); // avoid rounding error which can make E2 negative when M2 is negative return e2 > 0 ? e2 : 0; } /** Energy (timelike component of momentum-energy 4-vector) */ Scalar E() const { return std::sqrt(E2() ); } Scalar T() const { return E(); } /** vector magnitude squared (or mass squared) In case of negative mass (spacelike particles return negative values) */ Scalar M2() const { return ( fM >= 0 ) ? fM*fM : -fM*fM; } Scalar Mag2() const { return M2(); } /** transverse spatial component squared */ Scalar Pt2() const { return fPt*fPt;} Scalar Perp2() const { return Pt2(); } /** transverse mass squared */ Scalar Mt2() const { return M2() + fPt*fPt; } /** transverse mass - will be negative if Mt2() is negative */ Scalar Mt() const { Scalar mm = Mt2(); if (mm >= 0) { return std::sqrt(mm); } else { GenVector::Throw ("PtEtaPhiM4D::Mt() - Tachyonic:\n" " Pz^2 > E^2 so the transverse mass would be imaginary"); return -std::sqrt(-mm); } } /** transverse energy squared */ Scalar Et2() const { // a bit faster than et * et return 2. * E2()/ ( std::cosh(2 * fEta) + 1 ); } /** transverse energy */ Scalar Et() const { return E() / std::cosh(fEta); } private: inline static Scalar pi() { return M_PI; } inline void RestrictPhi() { if ( fPhi <= -pi() || fPhi > pi() ) fPhi = fPhi - std::floor( fPhi/(2*pi()) +.5 ) * 2*pi(); return; } // restrict the value of negative mass to avoid unphysical negative E2 values // M2 must be less than P2 for the tachionic particles - otherwise use positive values inline void RestrictNegMass() { if ( fM >=0 ) return; if ( P2() - fM*fM < 0 ) { GenVector::Throw ("PxPyPzM4D::unphysical value of mass, set to closest physical value"); fM = - P(); } return; } public: /** polar angle */ Scalar Theta() const { if (fPt > 0) return 2* std::atan( exp( - fEta ) ); if (fEta >= 0) return 0; return pi(); } // --------- Set Coordinates of this system --------------- /** set Pt value */ void SetPt( Scalar pt) { fPt = pt; } /** set eta value */ void SetEta( Scalar eta) { fEta = eta; } /** set phi value */ void SetPhi( Scalar phi) { fPhi = phi; RestrictPhi(); } /** set M value */ void SetM( Scalar mass) { fM = mass; if (fM <0) RestrictNegMass(); } /** set values using cartesian coordinate system */ void SetPxPyPzE(Scalar px, Scalar py, Scalar pz, Scalar e); // ------ Manipulations ------------- /** negate the 4-vector -- Note that the mass becomes negative */ void Negate( ) { fPhi = - fPhi; fEta = - fEta; fM = - fM; } /** Scale coordinate values by a scalar quantity a */ void Scale( Scalar a) { if (a < 0) { Negate(); a = -a; } fPt *= a; fM *= a; } /** Assignment from a generic coordinate system implementing Pt(), Eta(), Phi() and M() */ template PtEtaPhiM4D & operator = (const CoordSystem & c) { fPt = c.Pt(); fEta = c.Eta(); fPhi = c.Phi(); fM = c.M(); return *this; } /** Exact equality */ bool operator == (const PtEtaPhiM4D & rhs) const { return fPt == rhs.fPt && fEta == rhs.fEta && fPhi == rhs.fPhi && fM == rhs.fM; } bool operator != (const PtEtaPhiM4D & rhs) const {return !(operator==(rhs));} // ============= Compatibility secition ================== // The following make this coordinate system look enough like a CLHEP // vector that an assignment member template can work with either Scalar x() const { return X(); } Scalar y() const { return Y(); } Scalar z() const { return Z(); } Scalar t() const { return E(); } #if defined(__MAKECINT__) || defined(G__DICTIONARY) // ====== Set member functions for coordinates in other systems ======= void SetPx(Scalar px); void SetPy(Scalar py); void SetPz(Scalar pz); void SetE(Scalar t); #endif private: ScalarType fPt; ScalarType fEta; ScalarType fPhi; ScalarType fM; }; } // end namespace Math } // end namespace ROOT // move implementations here to avoid circle dependencies #include "Math/GenVector/PxPyPzE4D.h" namespace ROOT { namespace Math { template inline void PtEtaPhiM4D::SetPxPyPzE(Scalar px, Scalar py, Scalar pz, Scalar e) { *this = PxPyPzE4D (px, py, pz, e); } #if defined(__MAKECINT__) || defined(G__DICTIONARY) // ====== Set member functions for coordinates in other systems ======= template void PtEtaPhiM4D::SetPx(Scalar px) { GenVector_exception e("PtEtaPhiM4D::SetPx() is not supposed to be called"); throw e; PxPyPzE4D v(*this); v.SetPx(px); *this = PtEtaPhiM4D(v); } template void PtEtaPhiM4D::SetPy(Scalar py) { GenVector_exception e("PtEtaPhiM4D::SetPx() is not supposed to be called"); throw e; PxPyPzE4D v(*this); v.SetPy(py); *this = PtEtaPhiM4D(v); } template void PtEtaPhiM4D::SetPz(Scalar pz) { GenVector_exception e("PtEtaPhiM4D::SetPx() is not supposed to be called"); throw e; PxPyPzE4D v(*this); v.SetPz(pz); *this = PtEtaPhiM4D(v); } template void PtEtaPhiM4D::SetE(Scalar energy) { GenVector_exception e("PtEtaPhiM4D::SetE() is not supposed to be called"); throw e; PxPyPzE4D v(*this); v.SetE(energy); *this = PtEtaPhiM4D(v); } #endif // endif __MAKE__CINT || G__DICTIONARY } // end namespace Math } // end namespace ROOT #endif // ROOT_Math_GenVector_PtEtaPhiM4D